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1.  INTRODUCTION 

This  document  constitutes  the  final  report  for  Navy  Contract 
N00014-77-C-0583,  which  was  funded  by  NORDA  (Naval  Ocean  Research  and 
Development  Agency,  Code  500)  covering  the  period  from  11  July  1977  to 
1 January  1978.  The  object  of  the  research  performed  under  the  contract 
was  to  adapt  an  existing  radiation  transport  code  for  use  in  underwater 
acoustic  propagation  and  to  apply  it  to  a surface  duct  propagation 
problem.  The  new  computer  code  resulting  from  the  efforts  of  this 
contract  is  called  ACTOR  (Acoustic  Transport  in  an  Ocean  that  is  Random). 

During  the  past  three  decades  the  study  of  sound  propagation 

in  the  ocean  has  received  considerable  theoretical  and  experimental 

study.  The  theoretical  state-of-the-art  circa  1970  is  fairly  well 

summarized  in  the  books  by  Officer^  and  Urick^  which  also  contain 

references  to  a vast  amount  of  experimental  data.  Up  to  1970  the 

theoretical  studies  were  restricted  primarily  to  deterministic  oceans 

using  either  ray  trace  codes  or  normal  mode  codes.  This  does  not  mean 

that  fluctuation  phenomena  had  not  been  encountered.  As  long  ago  as 
(3) 

1947'  1 experimental  data  indicated  that  the  instantaneous  reverberation 

amplitude  is  Rayleigh  distributed.  To  treat  such  a problem  theoretically 

requires  detailed  oceanographic  data  as  well  as  advanced  numerical  methods. 

With  the  advent  of  the  parabolic  equation  (PE)  method  by  Tappert  and 
(41  (5) 

Hardin'  1 and  Tappert'  1 and  the  transport  equation  technique  by  Besieris 
and  Tappert, fluctuation  phenomena  can  be  treated  theoretically  in 
both  the  long  and  short  wavelength  regimes.  The  inadequacy  of  the 
oceanographic  data  is,  however,  still  a stumbling  block  for  high 
frequency  active  sonar  fluctuation  studies. 

In  the  low  frequency  regime  where  acoustic  propagation  over 
intercontinental  distances  can  be  considered,  the  characteristics  of  the 
deep  sound  channel  become  important.  It  is  now  generally  accepted  that 
fluctuation  phenomena  in  this  case  are  due  to  the  presence  of  internal 
waves.  For  many  regions  of  the  world's  oceans  the  space-time  scales  of 
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these  internal  waves  are  accurately  described  by  the  Garrett-Munk  model. 

This  model  of  internal  waves  in  conjunction  with  the  PE  method  has  been 
used  by  Flatte'  and  Tappert^^  to  describe  low  frequency  fluctuation 
phenomena  in  the  deep  sound  channel. 

At  higher  frequencies,  where  the  demands  of  a finer  space  grid 

make  calculations  using  the  PE  method  impractical  for  fluctuation  studies, 

the  transport  equation  finds  one  of  its  natural  areas  of  application. 

However,  because  the  transport  equation  is  completely  new  to  ocean  acoustics, 

we  felt  that  it  should  first  be  applied  to  an  area  where  experimental  data 

exist.  To  this  end  we  selected  the  surface  duct  propagation  study  reported 

by  Pedersen  and  Gordon. While  the  data  reported  in  Reference  11  goes 

out  to  10  kyd,  longer  range  data  to  about  65  kyd  was  generously  provided 
/ ip! 

by  Pedersen.  ' Of  particular  interest  to  us  was  the  cross  layer  propa- 
gation data  at  1030  Hz.  Beyond  6 kyd  neither  ray  theory  nor  normal  mode 
theory  can  adequately  account  for  the  level  of  ensonifi cation  observed  at 
the  400  ft,  below  duct,  receiver.  Using  the  transport  equation,  which  can 
easily  incorporate  scattering  by  rough  surfaces  and  by  internal  waves,  we 
investigated  the  possibility  that  these  energy  transfer  mechanisms  could 
account  for  the  energy  removed  from  the  surface  duct  and  detected  by  the 
receiver  below  the  duct. 

The  remainder  of  this  report  is  organized  as  follows.  The 
second  section  contains  a brief  discussion  of  the  transport  equation  and 
a rather  more  detailed  description  of  how  Monte  Carlo  methods  are  used  to 
obtain  its  solution.  The  various  scattering  kernels  used  are  taken  up  in 
Section  3.  The  experimental  data  and  a comparison  of  it  with  the  results 
from  normal  mode  theory  and  transport  theory  are  presented  in  Section  4. 

Various  conclusions  drawn  from  the  results  of  our  studies  are  given  in 
Section  5.  The  acknowledgements  of  the  persons  and  organizations  whose 
help  has  been  essential  for  the  successful  completion  of  this  work  are 
reserved  for  Section  6.  Finally,  the  references  are  provided  in  Section  7. 

Appendix  A contains  the  results  of  a preliminary  test  problem  in  which 
the  transport  code  with  no  scattering  mechanisms  was  compared  with  a ray 
theory  code  and  Appendix  B is  used  to  discuss  our  random  number  generator. 
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2.  THE  ACOUSTIC  TRANSPORT  EQUATION 


In  this  section  we  provide  a simplified  description  of  how  the 
transport  equation  can  be  obtained  from  the  parabolic  equation.  No 
attempt  is  made  at  mathematical  rigor  and  no  implication  is  implied  that 
this  is  the  only  way  to  obtain  the  transport  equation.  It  is  done  this 
way  only  to  show  that  there  is  a connection  between  the  transport 
equation  and  the  parabolic  equation  which  is  now  in  general  use  in  the 
underwater  acoustics  community. 


2.1  DERIVATION  OF  THE  TRANSPORT  EQUATION 

We  start  with  the  parabolic  equation  in  range  r and  depth  z 


i it  + 1 aft  + k 

ar  2k  ,22  2 


Ip  - 0 


Letting  0 = 0^  + <Sc$  we  can  put  Eq.  (2.1)  in  the  form 


(2.1) 


(2.2) 


where 


k = acoustic  wave  number  * ui/cQ 


ip  is  a complex  function  from  which  the  pressure  field  is  obtained  as 

P(r,z,t)  = f(r,z)e~1ut  . 

where 

f(r,z)  = <|>(r,z)H(kr) 
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» 


r 


r 

► > 


and  H(kr)  Is  the  outgoing  Hankel  function 

H(kr)~v^Fe1(kr"/4) 

appropriate  for  2-0  cylindrically  symmetric  solutions. 
Here 


V(r,z) 


(2.3) 


and  c^  and  cs  are  the  deterministic  and  stochastic  part  of  the  sound 
speed,  respectively,  and  cQ  is  a reference  sound  speed. 

This  is  mathematically  equivalent  to  the  Schroedinger  equation 
with  t -*■  r.  In  Eq.  (2.3)  6c$  is  a random  function  and  should  be 
written  6cs(r,z,a)  where  a given  o picks  out  a member  of  the  ensemble 
of  possible  <5cs  functions.  Now,  Eq.  (2.2)  can  be  written 

3*(r,z?J  , .2 

ii(/*(r,z1)  ^ [>/j*(r,z1)ip(r,z2)] 

az2 


+ V(r,z2)**(r,z1)*(r,z2)  = 0 . (2.4) 


4 


Also 

a**(r,z.)  , ,2 

1 1»  - H(r,z2)  — + ^ — 7 [<l»(r,z2)<|>*(r,Zj)] 

3Z1 

+ V(r,z1)\|;*(r,z1)K»(r,z2)  » 0 . (2.5) 


Subtracting  Eq.  (2.5)  from  (2.4)  and  defining 


(2.6) 


we  see  that  p satisfies  the  equation 

1 !r  + 5fH+  V(r*z2>a) 


For  propagation  primarily  In  the  r-dlrection  (l.e.,  small  e)  we  see  from 
Figure  2.1  that  the  wave  vector  in  the  z-direction  is  approximately  ke. 


k sine~ke 


. 
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If  we  now  introduce  the  Wigner  distribution 


(13) 


f(r,z,e)  • £ T dz'  e*-^02  <|>*(r,z+z' Mr,z-z' ) 

«/.ao 


(2.8) 


k ik0z' 

we  see  that  if  Equation  (2.7)  is  operated  on  by  - / dz'  e we  get 


i 2kez  * 


3( Z+Z 1 ) 3(Z -z'r 


+ £ /*  dz'  e12k0z  [v(r,z+z',o)  - V(r,z-z'  ,a)]  n»*(r,z+z‘  )n»(r,z-z' ) 

»-flD 

(2.9) 

Omitting  the  details  one  can,  by  an  integration  by  parts,  put  Eq.  (2.9) 
in  the  form 


31+  e 

3r 

3f 

3 A 

(to' 

2ikez 

» 

3 Z 

1 7T  J 

6 

k 

/■*■ 

2ik0z 

i TT  „ 

e 

(2.10) 


where 


kCd 


5c  (r.z.a) 


v a - k - 

v2  K c. 


(2.11) 


(2.12) 


The  function  f is  implicitly  understood  to  be  a function  of  the 
stochastic  part  of  the  index  of  refraction,  i.e.,  of  a.  Taylor 


r^r  r™ 


expanding  Vj(z  + z1)  about  the  point 
transport  equation 


|£+  8 
3r  3Z 


1 3 

7 37 


3f 

36 


we  can  obtain  the 


-j^/dz'  e2ik0z,[v2(z+z',a)  - V2(z-z\a)]**(z+z')*(z-z')  (2.13) 

where  the  entire  stochastic  aspect  of  the  problem  is  contained  in  the 
integral  on  the  right  hand  side.  Averaging  this  equation  over  the 
ensemble  of  possible  <5cs  Besieris  and  Tappert^’7^  show  that  the 
following  equation  results 


3F 

3r 


+ 


3JF 

30 


yd0'W(z,0,0‘  )[R(r,z,e' ) - F{r,zfe)]  - v(z)F 


(2.14) 


where  F(r,z,e)  * <f(r,z,0,a)>.  This  is  a transport  equation  for  small 
angle  scattering  where  the  scattering  is  determined  by  the  kernel 
W(z,0,e').  The  deterministic  part  of  tne  index  of  refraction  leads  to 
the  third  term  on  the  left  side  and  accounts  for  the  gradual  bending  of 
the  rays.  The  last  term  accounts  for  volume  absorption  with  depth 
dependent  absorption  coefficient  v(z)  (measured  in  inverse  range  units). 

The  details  of  the  scattering  kernel  W(e,e')  are  taken  up 
in  Section  3. 

The  physical  significance  of  the  function  F that  satisfies 
the  transport  equation  can  be  seen  as  follows.  From  Eq.  (2.8)  we  find 

I(r,z)  = J F(r,z,e)de  * < |<|>(r,z)  |2>  . (2.15) 


Thus  the  integral  of  F over  angle  gives  the  mean  acoustic  intensity  at 
depth  z and  range  r.  This  is  just  what  is  needed  to  compute  the  average 
transmission  loss.  Furthermore, 

J(r,z)  =y*F(r,z,e)ed0  = <-^■('1'*  ff  “ * |T)>  • (2-16) 

The  right  hand  side  is  just  the  mean  vertical  flux  of  acoustic  energy  as 
can  be  seen  by  writing  i|»  in  polar  coordinates,  ip  = Ae11*.  Defining  the 
geometrical  vertical  angle  as  the  vertical  phase  derivative,  9 = £ f|-  » 
we  obtain 

I(r,z)  = <A2(r,z)>  , (2.17a) 

J(r,z)  = <A2(r,z)e(r,z)>  . (2.17b) 

In  fact,  F contains  information  equivalent  to  the  mutual  coherence 
function  since  the  Fourier  transform  of  F with  respect  to  9 yields 

n(r,z,c)  s y*F(r,z,e)e'2ik9cde  = «p*(r,z  + ?)ip(r,z  - c)>  . (2.18) 

Therefore,  the  transport  equation  (2.14)  may  be  viewed  as  an  equation  for 
the  mutual  coherence  function.  As  the  range  r increases,  i.e.,  as  the 
separation  between  source  and  receiver  increases,  the  coherence  of  the 
acoustic  field  decreases  and  this  effect  is  just  what  Eq.  (2.14)  describes. 

Other  derivations  of  the  transport  equation  can  be  found  in 
References  14  - 17. 

2.2  VALIDITY  CONDITIONS  FOR  THE  TRANSPORT  EQUATION 

For  the  transport  equation,  as  given  by  Eq.  (2.14),  to  be  an 
accurate  representation  of  the  full  elliptic  wave  equation  in  a medium 
with  random  inhomogeneities  several  conditions  must  be  met.  These 


' 


"T^T 
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conditions,  as  summarized  by  Watson, are: 

1.  We  require  aX  « 1 where  a is  the  extinction 
coefficient  and  X is  the  acoustic  wavelength. 

Using  a wavepacket  picture  this  condition  says  that 
the  front  of  the  packet  should  not  undergo  significant 
extinction  relative  to  the  back,  i.e.,  the  packet 
should  not  suffer  severe  distortion  from  its  original 
form.  In  this  case  it  may  change  direction  but  still 
will  represent  a particle  with  given  momenta ; 

2.  The  eikonal  approximation  must  be  valid  in  the  absence 
of  internal  wave  scattering;  and 

3.  To  obtain  Eq.  (2.14)  the  long-time  Markovian  approxi- 
mation was  invoked.  This  requires  that  F(r,z,e) 
vary  slowly  over  the  horizontal  diffusion  length. 


2.3  THE  MONTE  CARLO  METHOD  OF  SOLUTION 

l» 

2.3.1  Introduction 

The  transport  equation  derived  in  Section  2.1  is  used  exten- 
sively in  many  areas  of  physics,  including  radiative  transfer  in  stellar 
atmospheres,  neutron  transport  in  nuclear  reactors,  and  to  calculate  the 
x-ray,  neutron,  and  gamma  outputs  from  nuclear  weapons.  Because  most  of 
these  applications  involve  complex  geometries,  realistic  cross  sections 
and/or  are  non-linear,  a vast  amount  of  effort  has  been  devoted  to 
developing  numerical  methods  for  solving  the  transport  equation.  In 
one-dimensional  geometries  the  predominance  of  effort  has  been  spent  on 
finite  difference  techniques,  many  of  which  are  described  by  Richtmyer 
and  Morton. 

For  two-  and  three-dimensional  problems  the  geometrical  com- 
plexities encountered  in  many  uses  has  led  to  the  development  of  sophisti- 
cated Monte  Carlo  methods.  As  applied  to  neutron  transport  the  Monte 
Carlo  method  is  aptly  described  in  Reference  19  while  the  method  as 
applied  to  non-linear  thermal  radiation  transport  is  discussed  by 
Fleck.'20'21’ 
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The  problem  of  calculating  acoustic  propagation  through  the 
ocean  is  complicated,  from  the  finite  difference  viewpoint,  because  of 
the  refractive  nature  of  the  deterministic  part  of  the  sound  speed  profile. 
This  is  compounded  by  the  fact  that  the  fluctuating  part  of  the  profile 
results  in  predominantly  forward  scattering  at  very  small  angles.  For  this 
reason  we  have  selected  the  Monte  Carlo  method  to  use  as  the  basis  for 
obtaining  numerical  solutions  of  the  acoustic  transport  equation  as  it  is 
represented  by  Eq.  (2.14).  It  is  easily  adaptable  to  complex  geometries 
and  complex  media.  In  addition,  the  Monte  Carlo  method  is  readily  modified 
to  accept  changes  or  additions  to  the  physics,  it  is  easy  to  incorporate 


time  dependence  into  it,  and  a multigroup  description  of  the  field  can 
easily  be  included. 

There  is  a vast  literature  available  concerning  the  many  facets 
of  the  Monte  Carlo  method  and  its  applications.  Subsections  2.2.2  and 
2.3.3  provide  only  a brief  glimpse  of  certain  important  features  of  the 
Monte  Carlo  technique.  Finally,  subsection  2.3.4  discusses  how  the  Monte 
Carlo  method  is  applied  to  the  specific  problem  of  acoustic  propagation  in 
a realistic  ocean  environment. 

2.3.2  Sampling  Techniques 

The  major  fundamental  difference  between  the  Monte  Carlo  method 
and  other  methods  of  obtaining  the  solution  to  a given  physical  or  math- 
ematical problem  Is  Its  reliance  on  random  sampling  to  represent  certain 
features  of  the  problem.  And  basic  to  all  random  sampling  techniques  is 
the  availability  of  a set  of  random  numbers,  £ , which  are  uniformly 

distributed  on  the  interval  (0,1].  The  random  number  generator  which  we 
use  Is  described  in  Appendix  B along  with  the  results  of  several  "goodness" 
tests. 

(a)  Function  Sampling 

To  construct  the  Monte  Carlo  solution  of  a transport  problem 
there  exists  a function  which  describes  the  probability  that  a Monte  Carlo 
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particle  will  go  a distance  x before  scattering.  The  probability  of 
the  particle  undergoing  a first  collision  between  x and  x + dx  is 
given  by 


p(x)dx  = e"aXodx  (2.19) 

where  o is  the  macroscopic  cross  section  of  the  medium  and  is  interpre 
ted  as  the  probability  per  unit  length  of  undergoing  a collision.  Inte- 
grating the  right  side  of  (2.19)  from  x = 0 to  x'  and  calling  the 
result  £ , gives 
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or, 


x' 


(2.20) 


By  sampling  £ uniformly  on  (0 , l] , Eq.  (2.20)  provides  the  correct 
probability  distribution  for  the  distance  to  collision.  This  example 
was  selected  because  it  is  used  later  in  sampling  the  distance  which  an 
"acoustic"  particle  or  ray  goes  before  scattering. 


In  general,  if  p(x)  is  a probability  density  function  on  the 
interval  a <.  x < b , then 

x 

C 3 P(x)  • f p(x)dx  (2.21) 

•'a 


determines  x uniquely  as  a function  of  £ . Often  an  analytic  solu' 
tlon  for  x cannot  be  obtained,  and  either  iterative  methods  must  be 
used  or  tables  of  partial  integrals  (obtained  numerically)  can  be 
employed. 
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(b)  Biasing  and  Importance  Samplin 

Biasing  and  importance  sampling  are  really  only  slightly  • 
different  aspects  of  the  same  thing.  A more  detailed  discussion  of 
these  ideas  can  be  found  in  Reference  22.  We  will  present  an  example  of 
biasing  and  importance  sampling  to  indicate  the  general  principles  in- 
volved. To  make  this  example  relevant  to  the  acoustics  problem  consider 
the  situation  shown  in  Figure  2.2.  Here  we  have  a region  of 


Source 

Z_ 


Region 

of 

Interest 


Figure  2.2.  Schematic  of  Biasing  in  Underwater  Acoustic  Propagation. 

interest  where,  say,  either  the  propagation  loss  and/or  the  angular  distri- 
bution of  the  acoustic  intensity  due  to  the  source  at  R * 0,  Z * l$  is 
reguired.  Since,  as  indicated,  only  rays  which  lie  in  the  two  bands  will 
pass  through  the  region  of  interest,  we  importance  sample  the  source  so 
that  the  vast  majority  of  rays  are  emitted  in  these  two  bands.  Then,  to 
account  for  the  overemphasis  on  the  number  of  rays  in  these  two  bands, 
we  bias  their  weights  to  properly  represent  the  energy  being  emitted  into 
these  two  small  angular  Intervals.  The  variances  of  the  answers  predicted 
in  the  region  of  Interest  typically  are  proportional  to  l/vR”  where  N 


«l 
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Is  the  number  of  rays  contributing  to  the  answer.  Hence,  by  increasing 
N we  reduce  the  variance. 

(c)  Estimators 

Typically  in  forward  analog  Monte  Carlo  four  types  of  estimators 
are  encountered,  generally  only  one  or  two  in  any  given  problem.  These 
are  called  the  collision  estimator,  the  last  event  estimator,  the  track- 
length  estimator,  and  the  point-detector  estimator.  The  acoustic  trans- 
port problem  has  two  cross  sections  which  characterize  the  medium.  These 
are  the  absorption  cross  section,  a , and  the  scattering  cross  section, 

a 

a , which  combine  to  form  a total  cross  section  oT  * + a . 

S I a S 

In  the  usual  acoustic  transport  calculation  the  ocean  is  divided 
up  into  zones  of  varying  sizes  (see  Figure  2.3)  which  all  have  the 
property  that  as  far  as  absorption  is  concerned  they  are  "optically"  thin, 
i.e.,  a ax  « 1 where  ax  is  the  maximum  dimension  of  the  zone.  To 

a 

calculate  the  acoustic  absorption  in  a given  zone  using  the  collision 


CM  CO 

II  II  II  • • 


Figure  2.3.  Typical  Zoning  for  an  Acoustics  Problem. 

estimator,  we  would  tally  E aa/ot  at  each  collision  which  takes  place  in 
that  zone  and  reduce  the  weight  E accordingly.  In  the  last  event  estimator 
the  particle  is  killed  when  e < oa/ot,  and  its  entire  weight  is  scored 
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in  the  zone  in  which  it  died.  If  the  random  number  c > °a/(7t  scattering 
occurs  and  no  weigh.t  is  lost.  The  collision  estimator  as  well  as  the  last 
event  estimator  lead  to  poor  statistical  results  in  optically  thin  regions. 
Because  the  acoustic  absorption  coefficient  is  so  small,  i.e.,  1 /oi  >>  ax 

u 

where  ax  is  the  maximum  zone  size,  we  must  consider  a better  estimator, 

the  path-length  estimator.  The  path-length  estimator  scores  the  quantity 
-oas 

E e in  the  zone  for  which  the  corresponding  path  length  is  s.  As 
discussed  in  Reference  22  the  relative  error  of  the  path-length  estimator 
is  much,  less  than  either  of  the  other  two  discussed  above. 

The  point  detector  estimator  would  not  be  particularly  useful 
because  of  the  refractive  nature  of  the  medium  and  hence  will  not  be 
discussed. 

2.3.3  Variance  Reduction 

(a)  Central  Limit  Theorem  and  Variance  Calculation 

The  variance  calculation  which  will  be  incorporated  into  the 

acoustic  transport  code  is  based  on  the  one  described  by  Cashwell  and 
(231 

I Everett'  ' which  makes  use  of  the  central  limit  theorem.  For  the  purpose 

of  describing  the  application  of  this  theorem,  let  us  c^afine  an  experiment 
E as  placing  a Monte  Carlo  particle  in  the  system  and  tracking  it  com- 
pletely out  of  the  system  (even  if  it  is  carrying  0 weight).  Now 
> define  a sample  space  ft  as  all  the  possible  histories  of  the  above 

initial  particle  placements  and  paths  through  the  system.  Also,  define 
a random  variable  on  this  sample  space  ft  as 

weight  of  particle  at  exit,  if  history  i of 
ft  exits  in  angular  group  a and  frequency 
group  v. 

0,  if  history  i of  ft  exits  in  some  angular 
group  a'  t a or  frequency  group  v'  t v. 

Assume  the  experiment  E is  performed  N different  times  (N  Monte 
Carlo  particle  trackings),  and  form  the  sum 


+ . . . + 


(2.22) 


XN 


Notice  that  the  are  mutually  independent  with  the  same  probability 
distribution.  If  we  let 


u * expectation,  or  mean,  of  (same  for  all  x^) 
2 

b ■ <3  * variance  of  (same  for  all  x^) 


(2.23) 


then  the  central  limit  theorem  states 


(ft-.l)  - -OS) 

(jSN  - N“|) 


as  N •*  * 


< Nue  = erf  / \ a$  N ^ » 

VM  / 


where  erf(x)  is  the  well  known  "error  function" 

Jl 


erf (x) 


-1  f e'c 
<7  Ja 


dt 


(2.24) 

(2.25) 


(2.26) 


which  is  tabulated  and  drawn  in  Figure  2.4. 

In  a radiation  transport  calculation,  we  are  interested  in  the 
probability  that  the  calculated  total  weight  exiting  the  system  in  a 
particular  angular  group  and  frequency  group,  as  determined  by  tracking  N 
Monte  Carlo  particles  through  the  system,  is  within  the  true  theoretical 
value  by  a factor  e.  We  can  use  the  central  limit  theorem  in  the  above 
form  to  give  us  an  estimate  of  this  measure  of  accuracy  by  making  the 
following  approximations  in  the  right-hand  side  of  Eq.  (2.25): 


1.  Replace  Nu  by  ; and 

2.  Replace  b by 


(2.27) 
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X ERF(X) 


0 

2227 

4284 

6039 

7421 

8427 

9103 

9523 


the  error  bars  to 

(2.28) 


A(dB  loss)  = 10  log1Qe  <un  0- 

= 4.3  . (2.29) 

Al  is  obtained  by  summing  over  the  product  of  the  particles  weight  and 
path  length  in  a zone.  Now  we  want  the  errors  to  imply  that  if  the  same 
calculation  (with  different  random  number  sequence)  were  repeated, the 
results  would  lie  within  the  error  bars  with  2/3  probability. 

So,  given 


In  the  acoustics  code  we  are  interested  in 
be  put  on  the  dB  loss  points.  Thus, 

dB  loss  = 10  log1Q 


and 


we  can  find  e such  that  erf^e/ZM-)  = 2/3.  With  this  e we  then  set 
the  error  bars  in  d8  loss  as  follows: 


A(dB  loss)  = 4.3(2e)  (2.30) 

Thus,  the  top  of  the  error  bar  is  at  dB  loss  - e(4.3)  and  bottom  is  at 
dB  loss  + e(4.3) . 

(b)  Splitting  and  Russian  Roulette 

The  usefulness  of  splitting  and  Russian  Roulette  is  very 
limited  in  acoustic  transport  applications.  The  primary  reason  is  that 
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the  typical  absorption  coefficients  for  frequencies  of  interest  are  so 
small  that  over  the  distances  of  interest  only  a small  fraction  of  the 
particles  weight  is  lost.  Furthermore,  the  distances  which  two  different 
particles  might  traverse  in  going  from  source  to  "detector"  do  not  differ 
significantly.  Thus,  all  particles  arriving  at  the  detector  do  so  with 
almost  the  same  weight.  Splitting  and  Russian  Roulette,  on  the  other 
hand,  are  most  useful  when  the  particles  arrive  with  a wide  spread  in 
weights.  In  such  a case  it  is  advantageous  to  kill  the  low  weight 
particles  or  particles  unlikely  to  contribute  to  the  answer  with  proba- 
bility p and  to  multiply  the  weight  if  it  survives  (probability  1 - p) 
by  the  factor  1/p.  In  a like  manner  a particle  of  weight,  E , which 
enters  a region  of  interest  might  be  split  into  n identical  particles 
each  of  weight  E/n. 

2.3.4  Simulating  Particle  Transport 
(a)  Source  Sampling 

The  solution  of  the  transport  equation  for  acoustic 
propagation  by  Monte  Carlo  methods  begins  by  sampling  the  source.  The 
sampling  procedure  is  shown  in  Figure  2.5.  We  assume  an  isotropic  source 


Figure  2.5.  Source  Sampling  and  Tracking  Schematic. 
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at  depth  Z$.  For  typical  problems  there  is  only  a limited  range  of  Ini- 
tial angles  for  which  energy  can  propagate  to  large  ranges.  This  range, 
as  denoted  above,  runs  from  +e  to  -e  relative  to  the  horizontal  where 
e * Arccos  (c$/cm).  Here  c$  is  the  sound  speed  at  the  source  and  cm 
is  the  maximum  sound  speed  along  the  profile.  Within  the  range  e to  -e 
rays  (particles)  are  started  equally  spacedin  angle.  Thus,  the  cosine  of 
the  Jtfl  particle  relative  to  the  vertical  axis,  C(J)  , is 

C(J)  3 cos  (ARG)  (2.31) 

where 

ARG  3 [(90  - 0 + irjjgg-  (2-32) 

end  N is  the  total  number  of  particles  considered.  Usually  3000  to 
5000  particles  are  used. 

(b)  Particle  Tracking 

As  a result  of  the  source  sampling  each  particle  begins 
its  trajectory  characterized  by  several  parameters.  These  are 

X 3 range  coordinate  = R 

Z = depth  coordinate 

A * cosine  of  trajectory  relative  to  horizontal 
C 3 cosine  of  trajectory  relative  to  vertical 
KI  3 K index  of  zone  in  which  the  particle  is  located 

LI  3 L index  of  zone  in  which  the  particle  is  located 

EN  3 weight  of  the  particle 

IE  3 flag  denoting  which  of  the  four  sides  of  a zone 
was  last  crossed. 

Each  zone  in  the  mesh  is  labeled  with  two  integers  K,L  as  shown  in 
Figure  2.3  and  A is  obtained  from  C as  A = . A is  assumed 

to  be  positive  since,  consistent  with  the  parabolic  equation  approximation 
we  allow  no  backward  reflected  energy. 
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Referring  to  Figure  2.5  we  have  a particle  entering  zone  0,  in 
which  the  gradient  of  the  index  of  refraction  is  (an/az)0  , at  point 
(r0.Z0).  Knowing  that  the  trajectory  in  zone  0 is  a circle  of  radius 
rQ  = l/(3n/3z)0  allows  a determination  of  the  center  of  curvature,  (R°»z°). 
to  be  made.  The  equation  of  the  circle  describing  the  trajectory  is  then 
solved  simultaneously  with  the  equation  for  each  zone  0 boundary.  After 
the  intersection  points  of  the  circle  with  each  of  the  four  sides  are 
known,  the  distances  from  (rq,Z0)  to  each  of  these  points  can  be  calculated. 
The  shortest  distance  to  boundary  (SDB)  is  then 

SDB  = MinfsJ1^  s[2),  S^,  S<2),  S<2),  sj!),  sj2)]  , 

where  sj1^,  s|2^  are  the  distances  to  the  two  points  where  the  circle 

crosses  the  ith  side.  If  any  si^,  is  less  than  10  ® it  is  set  to 

in 

10  . This  avoids  counting  the  point  where  the  particle  just  crossed  the 

boundary  (S^  = 0)  as  the  SDB.  After  the  particle  crosses  side  1 (see 
Figure  2.5)  its  new  direction  cosines  are  calculated  and  its  new  center  of 
curvature  (r*,Z*)  is  determined  using  the  new  value  of  (3n/3z)j.  The  new 
radius  of  curvature  is  r^  = l/(3n/3z)^.  This  procedure  continues  until 
the  particle  either  exits  the  mesh  at  the  maximum  range  of  interest  or  is 
assumed  to  be  absorbed  in  the  bottom  because  its  trajectory  crosses  a 
bottom  segment. 

(c)  Scattering 

In  the  process  of  tracking  the  particle  through  the  mesh 
it  could  either  scatter  from  the  surface  or  from  internal  waves.  If  a 
flat  surface  is  assumed  then  specular  reflection  obtains  and  at  the 
surface  the  C direction  cosine  merely  becomes  -C.  A rough  surface 
results  in  the  reflected  ray  coming  off  with  some  probability  distribution 
f (e i ,9 ) which  is  a function  of  the  incident  angle,  9.  , as  well  as  the 
exit  angle.  The  details  of  the  surface  scattering  probability  functions 
will  be  covered  in  the  next  section. 
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After  calculating  the  distances  to  the  various  zone  boundaries, 
as  discussed  above,  the  next  step  is  to  determine  the  distance  to  scat- 
tering. The  scattering  here  is  volume  scattering  and  internal  waves  are 
assumed  to  be  the  scattering  mechanism.  The  distance  to  scattering  is 
obtained  using  the  results  of  Eq.  (2.16), 

Ssc  s - ^ (2.33) 

where  o$  is  the  scattering  cross  section  and  £ is  a random  number. 

If 


the  particle  is  advanced  to  the  appropriate  boundary  as  discussed  before. 
On  the  other  hand,  if 


a scattering  takes  place  before  reaching  the  nearest  boundary.  In  this 
case  the  particle  is  advanced  along  its  trajectory  a distance  S to 
the  scattering  point.  Here  the  new  direction  cosines  are  determined  from 
the  scattering  probability  distribution  function  g(0,e'),  where  0'  is 
the  angle  relative  to  the  horizontal  after  scattering.  Given  the  new 
angle  0’  and  gradient  of  n , 3n/3z  , the  new  center  of  curvature, 

M M 

(R  ,Z  ) can  be  found.  Since  the  particle  has  not  changed  zones,  the 

L w 

radius  of  curvature  remains  unchanged.  With  these  new  conditions  an 
updated  trajectory  circle  Is  obtained  and  the  process  continues  as  before. 

The  details  of  the  internal  wave  scattering  kernel  will  also 
be  discussed  in  the  next  section. 

(d)  Transmission  Loss  Calculation 

To  calculate  transmission  loss  we  use  a method  based  on  the 
path-length  estimator  discussed  in  Section  2.3.2(c)  in  which  every 
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particle  which  passes  through  a zone  contributes  to  the  tally  in  that 
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zone. 


Consider  the  average  weight  of  the  itn  particle  in  zone  (K,L): 


_ SL  .•  -oa s' 

e^K.L)  = / eo  e ds'/cAt 
•/n 


(2.34) 


where 


tQ  = weight  of  the  particle  upon  entering  zone  (K,L) 
a = absorption  coefficient  i r.  zone  (K,L) 

a 

s^  = distance  travelled  in  zone  (K,L) 
c = velocity  of  sound  in  zone  (K,L) 

At  = time  step. 

The  quantities  c and  At  will  eventually  drop  out  of  our  expressions. 

If  we  assume  os.  « 1 , which  is  always  valid  for  the  a ‘s 

a i a 

and  zoning  sizes  which  we  use,  .then 


* -fsr 


(2.35) 


Averaging  this  over  the  volume  of  the  zone  gives 

.v  ^(K.L) 

ei(K,L)  = V0L(K,L) 


e1  s 
eo  5i 


(2.36) 


Here  and  R,,  are  the  two  R values  which  bound  the  zone,  aZ  * 

1^  - Zj  is  the  depth  thickness  of  the  zone  and  AR  = Rg  - R^  and  A$ 
is  the  azimuthal  angle  Increment. 
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Now,  at  a reference  range  of  r^  , select  a small  volume  of 

thickness  s ; then  the  reference  volume  is 
r 


Figure  2.6.  Schematic  of  Small  Volume  Element  at  Radius  rf. 
or 

2 

V = 2r  s sine  A<j>  . (2. 

r r r 

Hence 

- _ Esr  ,2 

eref  ~ cAt  Vr  1 

and  the  transmission  loss  becomes  (in  decibels) 


N(K.L)  • 10  logw  ^ 

r Zeis,  (2  slnejr?  I 


I 


where  the  sum  is  over  all  particles  which  pass  through  zone  K,L  and  E 
is  the  total  weight  emitted  through  the  reference  volume.  If  R and  Z 
are  measured  in  kilometers  and  the  dB  loss  is  reference  1 meter  then 
rr  * 10”3. 

The  final  step  is  to  assign  error  bars  to  the  values  of  N and 
this  is  done  as  discussed  in  Section  2.3.3. 


3.  SCATTERING  MOOELS 


In  this  section  we  discuss  the  various  scattering  models,  both 
surface  and  Internal  wave,  that  were  incorporated  into  the  ACTOR  code. 

For  surface  scattering  these  include  the  phase  screen  model  and  the  saw- 
tooth mode.  The  sinusoidal  model  was  studied  but  was  not  put  into  the 
code  because  it  appeared  to  differ  only  slightly  from  the  sawtooth  model 
at  least  as  far  as  the  average  scattering  angle  and  variance  are  concerned. 
This  will  be  discussed  more  fully  in  Section  3.1.3.  The  phase  screen 
model  is  developed  in  Section  3.1.1.  We  attempted  to  use  this  model  in 
conjunction  with  the  Pierson-Moskowitz^24^  spectrum  for  fully  developed 
seas;  however,  the  sea  state  which  obtained  during  the  experiment,  the 
results  of  which  we  were  attempting  to  calculate,  were  such  (extremely 
calm  sea)  that  it  did  not  apply.  Therefore,  we  proceeded  to  develop  the 
sawtooth  model  which  was  used  although  we  also  included  in  ACTOR  the 
coding  necessary  to  apply  the  phase  screen  model  should  it  be  needed  at 
some  future  time. 

The  internal  wave  scattering  models  are  discussed  in  Section  3.2. 
First  we  explain  how  the  results  of  the  Garrett-Munk^’^  model  are  incorpo' 
rated  into  the  ACTOR  code.  This  is  done  in  Sections  3.2.1  and  3.2.2.  We 
also  briefly  describe  another  model,  the  Gaussian  model,  in  Section  3.2.3. 

3.1  SCATTERING  BY  ROUGH  SURFACES 

3.1.1  The  Phase-Screen  Approximation 

As  pointed  out  by  Schnleder,^2^  the  treatment  of  scattering  by 
non-periodic  rough  boundaries  In  the  Monte  Carlo  method  cannot  be  done  by 
direct  simulation.  The  reason  Is  that  the  propagation  distance  scales  are 
so  large  compared  with  correlation  distance  associated  with  the  typical 
rough  ocean  surface.  Therefore,  we  simulate  instead  the  average  scattering 
"transfer  function"  which  takes  the  Incident  ray  with  angle  0Q  to  the 


exit  ray  with  angle  9 (see  Figure  3.1)  for  the  random  rough  surface 

described  by  the  random  function  c(r).  We  use  the  convention  that  upward 

going  rays  have  negative  angles  and  downward  rays  have  positive  angles; 

thus  < 0 and  9 > 0.  The  needed  relation  is 
o 


°r 

f(0,e,r)  = J G(0,eo)f(O,eo,r)d0c 

-CD 


(3.1) 


which  holds  for  0 < 0 < ® . 


For  a plane  wave  incident  at  angle  0Q  the  corresponding 


incident  pressure  is 


P.  = e 
tnc 


ik(r  cose„  - Z sine  ) 

O 0 


(3.2) 


Using  the  phase  screen  approximation,  near  Z = 0 the  scattered  wave  is 
given  by 


P . . e(k(r  cos9o  * 2 s1ne0>e1*(r) 

SC 


(3.3) 


The  total  pressure  field  is 


P = P.  + P 
me  sc 


(3.4) 


and  the  surface  boundary  condition  is 


P[r,5(r)]  = 0 


from  which  it  follows  that 


$(r)  = - 2kc(r)sine 


(3.5) 


The  scattered  field  at  Z = 0 therefore  is 


Psc(o,r)  = - e 


ik[r  cose  - 2e(r)  sine  ] 


(3.6) 
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In  the  parabolic  (small  angle)  approximation  we  write  this  as 


Psc(o,r)  = <j;(o,r)e 


(3.7) 


’l'(o.r)  = - e 


i$(r) 


(3.8) 


$(r)  = kr(cos0Q  - 1)  - 2ke(r)sin0Q  » - hkra^  " 2k£(r)0Q  . (3.9) 

Using  the  Wigner  function  to  calculate  the  transition  probability  G(e,eQ) 
we  get 

G(0,0O)  3 J*dy  e_l5lk9  y<«j>(r  - tyH*(r  + 4y)>  . (3.10) 


For  a surface  with  statistical  properties  given  by 


<52(r)>=  h2 


<S(r  + y)c(r  - >sy)>*  n p(y),  p(0)  » 1 


(3.11) 


and  assuming  c(r)  is  Gaussian  we  get 


k0rt  ~i k©_( 0 + ejy  -4k202h2[l  - p(y)] 


G(0’9o)  * 17  J dy  e 


o'  o‘J  a 0 
e 


(3.12) 


This  is  the  basic  formula  which  gives  G as  a function  of  the  correlation 
function  of  the  surface  fluctuations,  p(y).  Note  that 


yG(e,eo)d0 


1 . 


(3.13) 


Two  limits  can  now  be  investigated  based  on  the  Rayleigh  parameter, 

R =>  2khe  . * The  cases  R <<  1 and  R » 1 will  be  considered  and  then 
o 

combined  into  a single  Interpolation  equation. 


For  R » 1 we  use 


-4k292h2[l  - P(y)]  -2k^e^hV/L 
e as  e 


2„2h2..2/12 

o 


(3.14) 


where 


1 

7 


r! 

dy 


(3.15) 


y - 0 


The  quantity  L is  called  the  correlation  length  of  the  random  surface. 
The  rms  slope  is  given  by 


,2  *s 


<m>‘  ■ r 


(3.16) 


We  then  obtain 


-(e  + e )2L2/8h2 

S(0,0  ) - — e 
0 2h^T 


(3.17) 


The  rms  angle  of  scatter  is  found  to  be 

<(Ae)2>is  - ^ « 1 . 


(3.18) 


In  the  opposite  case  R <<  1 and  we  obtain,  skipping  some 


algebra. 


G(0,0O)  * (1  - 4k2eV)5(e  + e0)  + 0(0, 0Q) 


(3.19) 


where 


o <5  o ~7 * -iko_(o  9 )y 

2.2.2  o r . o o ^p(y)dy 


ke 


0(9, 9q)  * 4kt0jht  -£  Je 


(3.20) 


29 


The  first  term  on  the  right  in  Eq.  (3.19)  is  the  specular  term  and  the 
second  is  the  diffuse  scatter  term.  As  an  example  take 


2,,.  2 

p(y)  * /2L 


(3.21) 


Then 


, , , kL9  -k2L292(e  + 8 „)2/2 
D(e  + e ) « 4k2h2e2  — - e 0 0 

0 ° fa 


(3.22) 


The  width  of  the  diffuse  scatter  angular  spectrum  is  now 


20s  = 1 


<(A0)  = 


kL0, 


(3.23) 


Both  of  these  results  can  be  consolidated  into  a single  formula 
by  interpolation.  We  get 


G(0,0Q)  = S(0q)5(0  + 0Q)  + D(0,0q) 


(3.24) 


where  the  specular  component  is  given  by 


S(80)  * e 


-4k2h2e2 


(3.25) 


and  the  diffuse  component  is  given  by 

, -nr2^2 

" S 7 J — ? 


/ -4ktht0„\  , -(0  + 0„ 

0(0,0 J * yl  - e ) / 7 6 


)2/2a2 


(3.26) 


where 


2 1 . 4hfc 

0 ' TTT  T 
kVo^  l 


(3.27) 


30 


i 


r 


We  can  now  check  both  limits.  When  2kh9  >>  1 the  specular 

component  vanishes,  a = 4h<:/L  , and  D(e,eo)  reduces  to  the  previous 

expression.  On  the  other  hand,  when  2kh0  <<  1 the  specular  component 

2 2 2 2 0 

dominates  and  we  have  a = 1/k  L 9Q.  Again  we  recover  the  previous  result. 

Equation  (3.24)  is  sampled  statistically  in  such  a way  as  to  on 
the  average  reproduce  the  angular  distribution  implied  by  the  expressions 
for  S(9q)  and  D(9,0Q). 

To  understand  why  the  phase  screen  method  fails  when  used  in 
conjunction  with  the  Pierson-Moskowitz  spectrum  for  the  sea  state  (wind 
speed)  conditions  which  were  observed  during  the  experiment,  consider  the 
following. 

If  c(x)  is  the  wave  height,  let 


h2  =<52(x)> 


(3.28) 


then  from  Reference  24 


h2  = -SSH-  • 

4eg 

where  a and  B are  two  dimensionless  constants 


(3.29) 


a = 8.1  (10"'3)  , 

S * 0.74  , 

2 

u is  the  windspeed  and  g is  the  gravitational  constant  (9.75  m/sec  ). 
The  correlation  length,  L , is  then  given  by 


i - " 

<S*>H 


(3.30) 


where  the  average  slope  squared  is 


2 kc 

<S^>=  a Jin 


sg 


T 


(3.31) 
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Here  k is  a cut-off  wave  number  in  the  Pierson-Moskowitz  spectrum. 

c ( 26 ) 

One  way  to  determine  k£  is  to  consider  the  Cox  and  Munkv  ' expression 

for  the  average  slope  squared, 

S2  = (5.12  x 10"5  u ± .004) 


where  u is  the  windspeed  in  cm/sec.  At  u = 500  cm/sec  S':  = 0.0256  ± 


.004.  Taking  SQ  as  0.0256  we  get 


0.0256  = 8.1  (10‘3)  in 


k2  625 
c 


^0.74  (9.75)' 


or 


k = 1.63  m 
c 


-1 


So, 


<S2>  = 8.1  ( 10"3)£n(0.038  u4) 


(3.32) 


For  consistency  we  check  this  with  the  maximum  value  of  <S  > which  is 

taken  to  be  0.082  inferred  from  Kinsman^2'7^  who  claims  = 2/7  as 

the  theoretical  upper  limit  to  S.  Since  the  Pierson-Moskowitz  spectrum 

applies  to  fully  developed  seas  with  wind  speeds  of  20-40  knots  (10-20 

m/sec)  we  would  hope  that  when  u = 20  m/sec  is  substituted  into  Eq.  (3.32) 

2 

that  we  would  get  <S  > as  0.082.  Doing  this  gives 

max<S2>  = 8.1  (10“3)Jin  0.038  (20)4 
= 0.071 


which  is  reasonably  close  to  0.082. 


2 2 

h and  <S  > under  the  experi- 


Now  consider  the  results  for 
mental  conditions  of  u as  4 m/sec.  We  get  from  Eq.  (3.29) 
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P 

h2  = 0.0074  , 

and  <S2>  = 0.0184  from  Eq.  (3.32). 

For  the  surface  duct  in  which  the  measurements  were  made  the 
largest  angle  of  incidence,  0Q  , at  the  surface  was  about  3°  or  0.053 
radians.  The  Rayleigh  parameter  for  1000  Hz  sound  thus  becomes 

R = 2khe„ 
o 

3s  0.038 

where  k is  the  acoustic  wavenumber  34.2  m'1.  Since  the  Rayleigh 
parameter  is  very  small  we  have  almost  total  specular  reflection.  The 
real  trouble  arises,  however,  in  the  angular  distribution  of  the  diffuse 
component  the  width  of  which  is  given  by  Eq.  (3.23) 

<(48)2>^  . ^1- 

where 

L = h/<S2>'5  . (3.33) 

With  our  parameters  this  becomes 

<(a6)2>1s  * 7.13  radians  . 

Since  this  is  obviously  a nonphysical  result,  we  can  only  assume  that  the 
scattering  at  the  acoustic  frequencies  of  about  1000  Hz  or  less  is 
specular  and  takes  place  from  facets  which  result  from  the  swell  component. 
This  conclusion  is  what  led  us  to  develop  the  sawtooth  model  discussed 
next. 


r 


r 


3.1.2  The  Sawtooth  Model 

I Even  in  the  simple  sawtooth  model  the  actual  details  involved 

in  working  out  the  scattering  angle  for  all  cases  are  very  exten- 
sive. Therefore,  we  will  only  outline  the  method  and  provide  the  results 
for  the  various  situations  that  can  arise  when  an  acoustical  "packet"  or 
I ray  specularly  scatters  from  a sawtooth  surface  profile.  Thus  we  assume 

a surface  wave  of  the  form  shown  in  Figure  3.2, 


» 


[ 


> 


> 


I 


* 


* 


Incident  Packet 

Figure  3.2.  Schematic  of  Sawtooth  Surface  Wave. 


where 


S = 4aA 


(3.34) 


Geometrical  considerations  show  that  a ray  of  incident  slope  m1 
specularly  reflecting  from  a surface  with  slope  mg  will  exit  from  the 
reflection  with  slope  mg  given  by 


o 

2nVj  - mj  (1  - mg) 

* — 

1 + Zmgmg  - m^ 


(3.35) 


If  we  assume 


p 


iJ 


34 


i 


|m^|  « 1 ; | | « 1 . (3.36) 

* i.e.,  small  angles  of  incidence  and  shallow  waves,  then  Eq.  (3.35)  becomes 

s 2m.j  - m^  (3.37) 

t We  assume  the  packet  intersects  the  x-axis  at  xQ(0  <>  x < a),  and 

at  an  angle  9 (0  ^ 0 ^ir/2)  corresponding  to  a slope  m^  (0  ^ m^  < ®). 

We  need  to  distinguish  several  cases. 

Case  I,  x > A/2 
f 2 

This  is  the  case  shown  in  Figure  3.3. 


I Figure  3.3.  Geometry  of  the  Case  Where  xQ  > A/2. 

In  this  case  the  packet  always  hits  the  right  half  of  the  wave  (for  any 
positive  m^)  and,  with  m^  = - $ in  Eq.  (3.37)  , 

» 

m2  = - (2S  + m^)  . (3.38) 

This  equation  shows  that  the  slope  of  the  reflected  packet  is  always 
• negative,  no  matter  what  m^  is,  and  hence  the  scattered  packet  does 

not  interact  anymore  with  the  wave,  i.e.,  we  have  single  scattering.  It 


I 
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is  important  to  note  that  this  is  a consequence  of  assuming  small  angles 
(S  « 1 and  mj  « 1).  To  see  this  let  us  consider  the  general  reflection 
formula,  Eq.  (3.35).  If  the  incident  slope  were  given  by 

m.  = -»2S-  * (3.39) 

1 S*  - 1 

then  Eq.  (3.35)  predicts  m2  = 0 and  the  reflected  packet  would  travel 
horizontally  to  the  left  and  hit  the  left  side  of  the  wave,  resulting  in 
multiple  scattering.  But  for  S <<  1 , Eq.  (3.39)  has  no  solution  since 
we  restrict  to  be  positive  as  well  as  S. 

Case  II,  xQ  < \/2 

This  is  the  situation  shown  in  Figure  3.4. 


Here  there  are  two  subcases.  Depending  on  m^  , the  packet  may  hit  either 
the  right  half  of  the  wave  or  the  left  half.  Figure  3.4  shows  the  "cross- 
over" value  of  mj.  For  nij  less  than  that  shown,  the  packet  hits  the 
right  half;  for  m^  greater  than  that  shown,  it  hits  the  left  half. 
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Given  xQ  , let  us  compute  the  critical  value  of  m^  , which  we  call  m^. 
The  wave  surface  is 

y = Sx  (3.40) 


and  the  packet  trajectory  is 


y - n>i  (x  - xQ) 

Demanding  that  these  lines  intersect  at  x = X/2 


(3.41) 

leads  to  the  result 

(3.42) 


Case  II-A,  xQ  < X/2,  < m^ 

Here  the  packet  hits  the  right  half  of  the  wave  and,  from 
Eq.  (3.37), 


m2  = - (2S  + n^) 


(3.43) 


and  we  have  single  scattering. 

Case  II-B,  xQ  < X/2,  nij  > m^ 

In  this  case  the  packet  hits  the  left  half  of  the  wave  and 
Eq.  (3.37)  gives 


m2  = 2S  - m^  . (3.44) 

The  question  now  is:  Does  this  reflected  packet  hit  the  right  side  of 
the  wave  and  have  a second  scattering?  The  algebra  involved  in  answering 
this  question  is  lengthy  and  will  be  omitted.  The  result  is  essentially 
that  there  are  several  subcases. 
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Thus,  this  second-reflected  packet  cannot  go  back  and  hit  the  left  half 
of  the  wave  again.  Hence  a third  scattering  is  not  possible. 

Let  us  summarize  our  results.  Let  mi  be  the  slope  of  the 

incident  packet  (previously  called  mp,  let  mr  be  the  negative  of  the 

slope  of  the  scattered  (once  or  twice)  packet  and  S is  the  slope  of  the 

sawtooth  wave.  If  we  define  the  incident  angle 

e = tan~*  nu  (3.48) 

and  the  reflected  angle 

<t>  = tan'l  m (3.49) 

we  have  the  picture  shown  in  Figure  3.5. 


Figure  3.5.  Schematic  of  Incident  and  Final  Angles. 
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With  this  picture  we  then  have 


0 < XQ  < X/3 


mr  * mi  + 2S  » 0 <_  mi  < _ 2x  )S 


mr  ■ 4S  - mi  ■ (iTTr)s  i mi  i (rrr|s 

0 0 


mrsmi'2S  * (— ,Simi^ 

o 


X/3  < xQ  < X/2 


mr  = mi  + 2S  , 0 < mi  < (x  .-^X~)S 


mr ' "i  - 2S  • (r^r)s  i mi  * * 
• 0 


X/2  < xQ  < X 


m = m.  + 2S  , 0 < m.  < 

r i — i 


We  now  assume  that  the  phase  of  the  surface  wave  is  completely 

random  at  the  time  of  first  scattering  (and  that  the  wave  is  stationary 

until  the  packet  has  left  the  scattering  region).  We  want  to  ensemble 

average  our  results  over  all  possible  phases,  or  equivalently  over 

x„.  Thus 
0 


A 

"V = r j[K  "V 


(3.50) 


To  compute  the  variance  we  will  also  need 


fr  * f / dxo  "r 


(3.51) 


The  algebraic  details  of  these  calculations  will  be  omitted  and 
only  the  results  presented.  For  ny  we  have 


m.  < S 


rn  * m.  + 2S 
r 


(3.52) 


S <_  m.  j<  2S 


m * 4S  - — 
r m. 


(3.53) 


2S  <_  m.  <_  3S 


m 3 2m. 
r i 


6S  + 11  i- 
mi 


>_  3S 


(3.54) 


m„  * m.  + — — . 

r 1 m. 


(3.55) 


In  Table  3.1  we  have  listed  several  values  of  Z 3 mr/S  for 
various  values  of  5 3 m./S.  A corresponding  plot  is  given  in  Figure  3.6. 
“T 

For  mj;  we  get 


S 


m^  3 (m^  + 2S)^ 


(3.56) 
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Figure  3.6.  Mean  Value  of  Reflection  Tangent  mr  Versus  Incident  Tangent  mi  (S  = 4a/x). 


I 


» 


r 


t 


» 


» 


S < in.  _<  2S 


T ? 2 S 

4 - m,  - 2m,S  ♦ 16S^  - 6 ^ 

2S  < w.  <_  3S 


3 

17  = m2  + 2m.S  - 4S2  + 18 

>_  3S 

7 = m2  + 8S2 

r i 

2 

Defining  a variance,  o , as 

a2  s - (i.)2 

we  get 

<_  S 

O 

11 

CVJ 

O 

s < m1  i 

2S 

«2-"?(‘-s7)3(lts7) 


2S  <_  m.  <_  3S 

_3  <-4 

a2  = - 3m?  + 26m^  S - 84S  + 150  — - 121 


(3.57) 


(3.58) 


(3.59) 


(3.60) 


(3.61) 


(3.62) 


(3.63) 


i 


» 
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m^  >_  3S 


2 _ ..2  M S\ 
a - 4S  (1  - -j) 


(3.64) 


Note  from  this  equation  that  as  S goes  to  zero  so  does  a , as  it  should. 

2 2 

In  Table  3.2  we  tabulate  values  of  Z = a /S  for  various  values 
of  5 5 m../S  with  a corresponding  plot  given  in  Figure  3.7. 

The  results  just  obtained  are  used  in  the  ACTOR  code  by  fitting 
to  a Gaussian.  For  a given  incident  angle  e with  tangent  itu  we  obtain 
mr  by  sampling  from 


3.1.3 


P(^  ' 


The  Sinusoidal  Model 


("»r  - mr)2 


(3.65) 


For  the  swell  component  of  surface  roughness  we  used  the  saw- 
tooth model  in  our  calculations.  We  recognized  at  the  outset  that  this 
was  only  an  approximation  to  the  more  realistic  sinusoidal  character  of 
swell.  An  attempt  was  made  to  calculate  the  scattering  angle  as  a 

function  of  x and  m.  for  the  sinusoidal  model  as  was  done  for  the 
o 1 

sawtooth  model.  This  turned  out  to  be  analytically  impossible.  By  ig- 
noring second  scattering,  however,  we  were  able  to  obtain  analytic  so- 
lutions for  the  moments  of  mr  in  the  no-shadowing  region.  In  the 
shadowing  region,  again  ignoring  second  scattering,  these  moments  were 
obtained  numerically. 


The  no-shadowing  case  is  depicted  in  Figure  3.8.  This  is  the 


case  where 


mi  > 2,r  X 


(3.66) 


I 
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Figure  3.8.  Geometry  of  No-Shadowing  Case. 


and  if  we  use  the  same  definition  of  S as  in  the  sawtooth  model  (S  « 4a/ a) 
we  have 


m. 


(3.67) 


We  will  not  present  any  details  for  the  sinusoidal  model  but  will  only 
list  the  results.  For  the  no-shadowing  region  with  no  second  scattering 
the  first  moment  and  variance  become 


_2  ,2 

mr  = mi  + r ii^ 


2 

o 


2 <-2 
ir  S 

8 7 
mi  J 


(3.68) 


(3.69) 


The  case  of  the  onset  of  shadowing  is  shown  in  Figure  3.9.  For  xQ  values 

smaller  than  that  shown,  shadowing  becomes  important  and  results  for  the 

first  two  moments  of  m„  were  obtained  numerically. 

r 
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Incident 

Packet 

Figure  3.9.  Geometry  for  the  Case  of  Onset  of  Shadowing. 


The  sawtooth  profile  was  also  reinvestigated  assuming  no  second 
scattering  so  we  could  compare  with  the  sinusoidal  surface  in  the  non- 
shadowed  case,  Eqs.  (3.68)  and  (3.69).  The  results  are 


ny  (sawtooth)  = mi  + 


(3.70) 


ny  (sinusoidal ) 3 m. 


2 -2 
4 m,. 


(3.71) 


o2  (sawtooth)  = 4S2 


2 

a 


(sinusoidal ) 


(t  ) S2  1 


(3.72) 


(3.73) 


We  see  that  in  the  non-shadowed  region  the  sawtooth  model  is  brought  into 
agreement  with  the  sinusoidal  model  if,  in  the  sinusoidal  model,  we  set 
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(3.74) 


FT  " . . ‘ ' ‘ ' J y IUII 

. 


I 


In  Figure  3.10  we  show  ny/S  for  the  sawtooth  with  no  second 
scattering  plotted  against  m^/S.  While  in  Figure  3.11  we  have  the 
corresponding  results  for  the  sinusoidal  model.  Note  that  in  this  case 
we  plot  mr/S’  versus  m./S*  where 


» 


2tt  a 


whereas  for  the  sawtooth  model  S is  used,  where 


(3.75) 


S = 4 ~ . (3.76) 

This  transformation  makes  the  model  results  identical  in  the  non-shadowed 

2 

region.  Figures  3.12  and  3.13  contain  the  corresponding  results  for  a . 


9 


9 


9 


First,  we  compare  the  results  of  the  sawtooth  model  with  and 
without  second  scattering,  i.e.,  Figures  3.6  and  3.10.  We  see  the 
following  effects  due  to  the  neglect  of  second  scattering: 

1.  For  m^  <_  S there  is  no  effect  on  either  the 
mean  or  the  variance; 

2.  For  m^  >_  3S  there  is  no  effect  on  either  the 
mean  or  the  variance; 

3.  For  S _<  m^  <_  3S  the  neglect  of  second  scattering 

(a)  decreases  the  mean 

(b)  increases  the  variance. 


The  surface  duct  problem  that  will  be  discussed  in  Section  4 had 
S * 0.0025  with  mi  ^ tan  2°  * 0.035  so  m^/S  % 14.  Thus,  in  the  sawtooth 
model  for  this  problem  we  can  safely  neglect  second  scattering. 
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Figure  3.10.  Mean  Value  of  Reflection  Tangent  mr  Versus  Incident  Tangent 
m.| , Neglecting  Second  Scattering  for  Sawtooth  Model 
(S  =■  4a/ x) . 


We  now  compare  the  sinusoidal  model  without  second  scattering 

with  the  sawtooth  model  also  without  second  scattering.  In  the  non- 
_ 2 

shadowing  region  mr  and  a as  a function  of  S'  have  the  same  form 
as  in  the  sawtooth  model  when  plotted  as  a function  of  S.  Thus,  for 

m.  > 2*  7 = /I  S' 

1 A 

the  sinusoidal  and  sawtooth  models  are  in  exact  agreement  if  S'  in  the 
sinusoidal  model  replaces  S in  the  sawtooth  model.  For  the  case 

m^  < S' 

there  is  a smoothing  which  occurs  in  the  sinusoidal  model  due  to  the  fact 
that  a distribution  of  scattering  slopes  is  available.  Also,  in  the 
sinusoidal  model  as  m./S'  approaches  zero  so  does  mr/S‘.  This  is 
because  the  available  points  on  the  sinusoid  where  scattering  can  take 
place  are  those  corresponding  to  very  small  slopes. 

We  conclude  that  for  our  applications  the  sawtooth  model  is 
perfectly  adequate.  If,  in  the  future,  situations  are  encountered  where 
m./S'  is  very  small,  provisions  would  have  to  be  made  to  more  accurately 
include  the  region  below  m./S'  3 . Also,  the  effects  of  multiple 

scattering  would  have  to  be  investigated. 

3.2  SCATTERING  BY  INTERNAL  WAVES 

3.2.1  The  Total  Scattering  Cross  Section 

In  this  section  we  discuss  how  the  internal  wave  spectrum  of 
Garrett  and  Munk^®*^  is  sampled  for  use  in  the  Monte  Carlo  acoustic 
transport  code,  ACTOR. 

From  Reference  28  the  acoustic  scattering  cross  section  can 

be  written 

4 

a(q,q')  = ft  F(k,Q,Z)  (3.77) 


55 


r 


» 


» 


i 


» 
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where  F(k,Q,Z)  is  the  power  spectrum  of  sound  speed  fluctuations  and 
the  following  definitions  hold: 

A 

q'  = unit  vector  in  direction  of  incident  packet 

A 

q * unit  vector  in  direction  of  scattered  packet 

K = q(q  - q* ) 

q = wavenumber  corresponding  to  momentum  transfer 

vector,  K 

K = (k,Q) 

n 2~  -*■ 

k = /kj^  + k2  = horizontal  component  of  K 
Q = vertical  component  of  K 

We  are  interested  in  evaluating  the  scattering  coefficient 

a$(q’)  =Ja(q,q,)d(iq  . (3.78) 

Since  we  are  integrating  over  all  q directions  with  q'  fixed  we  see 


Figure  3.14. 


JSeometry  Schematic  for  dnq 


Calculation. 
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from  Figure  3.14  that 


dflq  = cosy  dyc4 


dyd<>  for  small  y . i.e.,  small  vertical  angles.  (3.79) 


dK  - (dk.dQ) 

= ( qd4>  ,qdY ) 


(3.80) 


q 


(3.81) 


.(q* ) - y§^  F(k.Q,Z)  ^ 


= 4/f  F(k,Q,z)dkdQ  0 < k < Q . 

Using  the  Garrett-Munk^8,9^  power  spectrum  we  get 

SW)  ■ fq2  T dQ f dk  (3.82) 

min 


where 


8 2 <4*i 

a “ 7 <v  > ^ r 


2 “M 

SZ  * -J  K<  1 
n 


(3.83) 


(3.84) 
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Q • s a minimum  allowable  change  in  vertical  wavenumber 
upon  scattering  due  to  finite  depth  of  the  ocean 

Q*  * J**  n/(BnQ)  (3.85) 

Q = j it  n/(BnQ)  . (3.86) 

In  the  above  equations  the  following  definitions  hold: 


<u2>  * <Wg>  (n/nQ)3 
2 

<u  > * variance  of  the  sound-speed  fluctuations  extrapolated 
0 to  the  surface 

5 2.5  x 10'7 

j 3 inertial  frequency  (7.3  x 10’5  sec"1  at  30°  latitude) 

n 3 n(Z)  3 buoyancy  (Brunt-Vaisala)  frequency 
B 3 scale  of  the  sound  channel  (~  1 km) 

-3  -1 

nQ  3 n at  the  surface  (5.2  x 10  sec  ) 

j*  3 constant,  suggested  in  Reference  29  to  be  3. 

Equation  (3.82)  can  be  written 


os(q') 


fq2  f dQ 

%1n 


vs- 

(Q2  + Q2) 


/ 


dk 


[7T7W 


(3.87) 


Letting 


* 


we  get 


k = 0 =£>  x = e2Q2 


k = Q =>  x = (1  + e2)Q2 


as(q')  55  f q2  T dQ  — 2“^ 

s 4 o/  (0^  + 


Qfon  (Q  + Q*)  6 Q^(  1 + S )Q" 


and  using  the  fact  that  e « 1 we  get 


oA q')  - 


f — 

48  %n  ^ + «> 


(3.88) 


(3.89) 


We  must  now  assign  a value  to  Qmin.  Presumably  2w/(ocean  depth )<_  Q*. 
Slnce  ^min  only  aPPears  in  the  argument  of  a log,  its  actual  value  is 
rather  unimportant  and  we  choose  Qmi-n  - Q*.  We  then  get 

«.($'>  — ZjXmd) 
s 4r  2q;  1 


= -4  q2  <u2>  — A-  in(2)  . 


J* 


(3.90) 


Using  j*  = 3 , B = 1.2  km  and  scaling  to  a depth  of  -B  gives 


* 013  42  & 


'•56  4z 


dB/km 


(3.91) 
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where  f is  the  acoustic  frequency  in  kilohertz  and  the  scaling  used  the 
canonical  form 


n(Z)  - noeZ/B,  Z < 0 


(3.92) 


3.2.2  Sampling  the  Scattering  Kernel 

After  integrating  over  the  horizontal  wavenumber,  i.e..  in  going 
from  Eq.  (3.87)  to  (3.89),  we  obtain  the  change  in  the  vertical  component 
of  K - q(q  - q')  to  be  given  by  the  density  function 


in  which 


p(Q) 


1 

Q(Q2  +~  Q*) 


p(Q)dQ 


-V 

2Q ; 


(3.93) 


(3.94) 


To  sample  from  p(Q)  we  normalize  by  the  right  side  of  Eq.  (3.94)  and 
using  the  results  of  Section  2.3.2  we  set 


rQ  dQ1 

Q'(Q'2  + Qt) 

\ tn(2) 

2q; 


(3.95) 


or 


C tn(2) 


(3.96) 
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•ft 


ft 


» 


I 


I 


> 


hence 

Q* 

Q - + - — ■ (3.97) 

J 2C  - 1 

When  c * 1,  Q » ± Q*  as  it  should  and  in  Table  3.3  we  list  various 
values  of  Q for  the  corresponding  $. 


Table  3.3.  Values  of  Q Corresponding  to  Various  Values 
of  the  Random  Number,  £. 


5 

Q 

1.0 

± Q* 

0.75 

± 1.21  Q. 

0.50 

± 1.55  q* 

0.25 

± 2.30  q* 

0.10 

± 3.73  q* 

0.01 

± 12.00  q* 

0.001 

± 38.00  q* 

from  this  table  that  | 

q|  satisfies  the 

fol lowing 

50$  of  the  time  | Q | 

will 

1 Je  between 

Q* 

and 

1.55  q* 

75$  of  the  time  | Q | 

will 

lie  between 

Q* 

and 

2.30  q* 

90$  of  the  time  ( Q ( 

will 

lie  between 

Q* 

and 

3.73  q* 

99$  of  the  time  | Q | 

will 

1 1e  between 

Q* 

and 

12.00  q* 

99*9$  of  the  time  |q| 

will 

lie  between 

Q* 

and 

38.00  q* 

Translating  these  Q's  Into  changes  In  vertical  angle  we  get 
(at  Z ■ -B  ■ - 1.2  km) 

Q*  - - 2.88  km-1  . (3.98) 
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Now,  the  acoustic  wavenumber  is 


q = — = 4 X 103  km"1  at  f*  1000  Hz 
c 


so,  referring  to  Figure  3.15  we  see  that 


Figure  3.15.  Vertical  Angle  Change  Due  to  Scattering. 


for 


Q - Q* 

Q = 4 Q* 

Q - 12  0* 


A6  = 0.04 
A0  = 0. 16c 


A0  = 0.48 


We  assume  small  angles  and  that  A(Tane)  = A0. 

An  approximate  picture  of  the  probability  density  function  is 
shown  in  Figure  3.16. 

Watson^®)  derived  an  expression  for  the  n ns  vertical  scatter 
angle  (in  radians)  for  a ray  on  the  sound  channel  axis.  He  found 


“T  *5 


(0")  8 1.1  x 10-3  R^m 


(3.99) 
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The  derivation  made  use  of  the  diffusion  approximation  in  angle  in  which 
no  limit  is  placed  on  the  smallness  of  je  - e'|. 

We  simulated  the  same  problem,  i.e.,  a channel  axis  ray  under- 
going multiple  scattering  as  it  proceeds  in  range,  using  Monte  Carlo 
methods  in  which  the  scattering  cross  section  was  obtained  from  Eq.  (3.91) 
and  the  scattering  angle  was  obtained  by  sampling  from  Eq.  (3.98),  the 
probability  density  of  which  is  shown  in  Figure  3.16.  An  interesting 
difference  between  the  Monte  Carlo  method  and  the  diffusion  method  lies 
in  the  fact  that  the  Monte  Carlo  method  samples  from  a distribution  which 
requires  each  scattering  angle  to  exceed  some  minimum, while  the  diffusion 
method  allows  all  | e - 0 ' | >0.  It  is  remarkable,  however,  that  both 
methods  give  essentially  the  same  results  for  the  rms  scattering  angle 
as  a function  of  range  as  can  be  seen  from  Figure  3.17.  This  might  not 
be  surprising  at  large  ranges  after  many  scatterings  have  occurred,  but 
the  good  agreement  even  at  10  km  (~  2 scatterings)  is  surprising. 

We  did  not  check  the  angular  distribution  as  a function  of 
range  but  it  would  not  be  surprising  if  substantial  differences  appeared 
between  the  results  predicted  by  the  two  methods  at  short  ranges.  Some 
effort  will  be  made  to  investigate  this  at  a later  date. 

3.2.3  The  Gaussian  Model  for  Internal  Wave  Scattering 

In  Sections  3.2.1  and  3.2.2  the  Garrett-Munk  model  for  internal 
wave  scattering  was  discussed.  A salient  feature  of  that  model  as  pre- 
sented was  the  neglect  of  range-dependent  sound-speed  fluctuation  corre- 
lations. This  led  to  a scattering  kernel  which  depended  only  on  the 
difference  between  the  incident  and  scattered  angles  and  the  scattering 
Integral  becomes  a convolution  Integral. ^ The  development  by  Watson'2®^ 
also  ignored  range-dependent  correlations.  We  will  show  here  how  the 
Garrett-Munk  formalism  changes  when  a more  accurate  expression  is  used 
for  the  horizontal  change  in  wavenumber.  This  results  from  taking  into 
account  the  range-dependent  sound-speed  fluctuations.  We  will  then 
suggest  a simpler  model,  the  Gaussian  model,  for  doing  approximately  the 
same  thing. 


00 


Figure  3. 


Range  (km) 


RMS  Vertical  Scattering  Angle  for  a Channel  Axis  Ray  as 
Determined  by  Diffusion  in  Angle  and  by  Monte  Carlo. 


Using  9 as  the  small  vertical  scattered  angle  with  e the 
incident  angle  we  have 


*•■>■//£  F(k,e,z,)ded<j> 


Q = q(e-e') 


\ = U-<t>')2  + (He2-^'2)2 
q 


(3.100 


(3.101 


where  the  notation  is  the  same  as  in  Section  3.2.1.  It  is  the  second  term 

2 

in  this  expression  for  ( k/q ) that  is  new.  Thus, 

o$  g | 1 1 fr(q4<t>-<i>'  )2+(}502-j50|2)2.  q(e-0'),z)ded<|> 

» J" g(e,9‘,z)de  (3.102 


where 


g(9,9',z)  = tWu-*1)2  + (^e2-^'2)2,  q(e-e'),z)d* 


q | 9-9  1 


2!  oq 


3 1 0-9' 


[q2(0-9')2+Q2*]  JQ  [k2+q2(He2-Js0 ’ 2)2+62q2(0-9 ' ) 


After  performing  the  integration 


9(0,0' ,z) 


|8-8'|[(e-e')V.(z)]  [s2(z)  ♦ (^-)  ] . 


where 


Q*(z)  ■ (J*n/qB)  (n(z)/n  ) 


and 

g(z)  = o»i/n(z) 

The  function  g is  the  rate  per  unit  distance  along  a ray  for  scattering 
from  6'  to  0 at  depth  z.  It  is  even  in  0 and  0 1 , which  is  the 


1 1 


principle  of  detailed  balance. 

The  additional  term,  not  found  before,  is  [(e+e* )/2]  . Note 
that  it  is  only  important  when  (0+0* )/2  > s(z).  The  parameter  0 measures 
the  degree  of  anisotropy.  If  0 were  of  order  unity,  then  the  additional 
term  would  not  be  needed  because  all  angles  are  small  anyway.  The  fact 
that  corrections  may  be  needed  for  angles  of  order  of  the  anisotropy 
factor  is  here  confirmed  explicitly.  We  have 

0 ( z ) = ■—  e‘z/B  , z < 0 . 

0 

For  the  canonical  values,  0(z)  3 0.014  e"z  where  z is  in  kilometers 
and  Table  3.4  lists  the  values  of  0 and  0*  for  various  depths. 

Table  3.4.  Depth  Dependence  of  0 and  0* 


z(km) 

0(rad) 

0(deg) 

0*(rad) 

e*(deg) 

0 

0.014 

0.8 

0.00225 

0.13 

-1 

0.038 

2.2 

0.00083 

0.05 

-2 

0.104 

5.9 

-3 

0.282 

16.2 

Also,  at  1 kHz, 

8*(z)  3 

0.00225 

ez  . 

and  it  is  seen  that  the 

correction 

term  is 

only  needed  near  the 

where  the  anisotropy  is  large. 
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The  total  cross  section  is  given  by 


r(e'.z)  - I 


o r 09 

J | e-e ' | [(0-0')2+0*2  [e2+(^-)  _ 


0'-0* 


- ™ - ■ ■ - 2“ 

|0'-0|  [(0-0')2+0*2][e2+(^-)  ] . (3.105) 


This  can  be  put  into  the  form 


oo 

cr(e ' ,z)  3 §■  / dsj|-  * "f — ) f " J + " j}  ' 

V r+l)  (s2+(e.4*5.)  82+( e ■ ) 

which  is  much  more  complicated  than  before. 

In  the  Gaussion  mode 1 the  scattering  kernel  is  given  by 


w(z,0,0‘)  = 2Trq?(z,0-0'  , ^O2-^'2) 


(3.106) 


(3.107) 


where 


fU.e.u)  > — /dz'  f dr'e“lq^ez,“ur  ^r(z,z‘ ,r) 

J J 


(3.108) 


r(ztz',r')  * -5-<6c(z+J5Z,,r-Hsr,)><5c(z-J5z' ,r-igr')>  . (3.109) 


To  say  more  about  W,  one  needs  to  know  the  correlation  function 

A 

of  sound  speed  fluctuations  r(z,z',r')  or  its  power  spectrum  r(z,0,u). 
Usually  one  is  only  able  to  estimate  the  horizontal  and  vertical  scale 


1 


f 


lengths  (LH  and  Lv,  respectively)  and  the  depth  dependence  of  r(z,o,o)  * 
1/Cg  <<5c(z,r)<5c(z,r)>.  It  is  generally  believed  that  r(z,o,o)  is  largest 
near  the  main  thermocline  and  decreases  rapidly  below.  Thus  a reasonable 
but  crude  model  for  r(z,z',r')  which  contains  only  four  parameters  is 


r(z,z‘  ,r*)  = 


2 "Z/BQ 
a e e 


-z 


,2 


/2LJ  -r 


,2 


/2L 


-4 

Typical  values  might  be  o»5xl0  , B=*l  km,  Ly^lOO  m,  and  LH=*10  km. 
An  expression  for  W then  follows  immediately 


W( Z , 8 , 9 1 ) 


L L k3  2.-3Z/2B  '^o'9-9')2  ■*&*«*) 
LV  H o e e e 


where  a = (4c/cQ). 

This  is  the  Gaussian  model  alluded  to  earlier.  It  contains, 

in  at  least  an  approximate  manner,  the  essential  oceanography  and  we  feel 

it  would  be  very  useful  in  performing  internal  wave  scattering  sensitivity 

studies.  We  have  not  as  yet  incorporated  it  into  ACTOR;  however,-  at  a 

later  date  when  comparisons  of  transport  theory  with  the  method  of  path 
(291 

integrals'  ' are  made,  we  expect  the  Gaussian  model  to  play  an  important 
role  precisely  because  it  does  account  for  range-dependent  sound-speed 
fluctuation  correlations  whereas  the  method  of  path-integrals,  as  applied 
to  date,  does  not. 
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4.1 


4.  THE  EXPERIMENTAL  DATA  AND  ITS  COMPARISON  WITH 
TRANSPORT  AND  NORMAL  MODE  CALCULATIONS 


THE  EXPERIMENTAL  ACOUSTIC  DATA 

The  experimental  data  and  the  normal  mode  calculations  with 


which  we  will  compare  our  transport  results  were  both  supplied  to  us  by 

f og  \ 

Mel  Pedersen  and  DeWayne  White'  ' of  the  Naval  Ocean  Systems  Center, 

San  Diego,  California.  They  also  provided  the  expressions  for  the  sound 
speed  as  a function  of  depth. 

The  following  description  of  how  the  experimental  data  were 
obtained  was  taken  from  Reference  31. 

"The  experimental  data  were  obtained  in  1955  off  the 
California  coast  in  2200-fathom  water.  One  ship, 
towing  two  acoustic  sources,  opened  range  from  a 
second  ship  on  a fixed  course  at  a speed  of  about 
three  knots.  These  two  sources  were  530-  and  1030-cps 
Fessenden  oscillators,  which  were  strapped  together 
rigidly  and  towed  at  a depth  of  about  54  ft.  This 
depth  was  monitored  on  a pressure  (depth)  indicator. 

The  two  sources  were  pulsed  alternately  with  evenly- 
spaced  500-msec  pulses.  Eight  pulses  per  minute 
were  transmitted  at  each  frequency. 

Signals  were  monitored  on  hydrophones  at  nominal 
depths  of  50,  232,  and  400  ft.  The  sound  pulses  at 
the  two  frequencies  were  separated  by  means  of 
filters  and  were  recorded  on  a six-channel  recorder. 

Acoustic  travel  times  were  obtained  by  measuring  the 
time  difference  between  a pulse  received  by  radio 
transmission  and  the  corresponding  acoustic  pulse 
received  on  the  hydrophone  at  50-ft  depth." 

From  this  description  we  see  that  one  pulse  was  received 
during  each  10-yard  range  interval.  The  first  step  in  reducing  the  vast 
amount  of  data  was  to  average  it  over  100-yard  intervals.  Thus 


t > 


V!ref 


(4.1) 


► » 
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(d)2^I(v. 


(4.2) 


The  transmission  loss  curves  are  then  obtained  by  plotting  10  log^ 
(T/Iref)  versus  range. 

The  transport  equation  results  provide  an  expected  value  for 
I / 1 re ^ averaged  over  all  members  of  the  internal  wave  ensemble.  The 
measurement,  on  the  other  hand,  basically  picks  out  only  one  member  of 
the  ensemble  since  only  one  measurement  was  made  at  each  ten-yard  interval 
and  the  time  between  measurements  was  small  compared  with  an  internal 
wave  period.  If  we  assume  that  internal  wave  phenomena  represent  a 
stationary  random  process,  then  we  can  look  at  the  time  series  in  range 
to  deduce  a range  averaged  variance. 


■■  (d!  ■ te)! 


(4.3) 


To  be  sure,  the  actual  standard  deviation,  a , is  a function  of  range 
being  smaller  in  regions  of  high  intensity  and  larger  in  regions  of 
low  intensity. 

(32 

The  standard  deviation  data  were  provided  to  us  by  Mel  Pedersen' 
as  a function  of  range  for  each  frequency  and  for  each  depth.  We  will 
not  present  that  data  but  only  provide  average  values  of  a "^'2  for 
each  case.  These  average  values,  after  appropriate  scaling,  are  given 
in  Table  4.1  and  are  to  be  interpreted  as  meaning  that  the  vast  majority  of 
the  range  dependent  o's  for  each  range-depth  case  lie  below  the  corres- 
ponding value  listed  in  the  table. 
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Table  4.1.  Average  Standard  Deviation  in  Decibels 
as  a Function  of  Range  and  Frequency. 


Frequency 

(Hz) 

Depth 

(ft.) 

a 

(db) 

530 

50 

0.7 

530 

232 

0.8 

530 

400 

2.0 

1030 

50 

0.8 

1030 

232 

2.5 

1030 

400 

4.0 

4.2  THE  SURFACE  DUCT  MODEL 

The  vast  majority  of  the  transport  and  normal  mode  results  pre- 
sented in  the  next  section  were  obtained  using  the  surface  duct  model 
shown  in  Figure  4.1.  This  is  a two-layer  model  with  layer  number  one 
being  the  top  layer.  The  index  of  refraction  in  the  ith  layer  is  given 
by 

n*  » 1 - 2 ^ (Z  - Z.)  (4.4) 

where  Z.  is  the  depth  at  the  top  of  the  ith  layer,  c,  is  the  sound  speed 
at  the  top  of  the  i u layer,  and  characterizes  the  strength  of  the 
gradient  in  the  i layer.  The  parameters  for  the  two-layer  model  are 
given  in  Table  4.2. 


i 
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Table  4.2.  Parameters  Characterizing  the  Two-Layer  Model. 


Layer 
Number,  i 

ci (kyd/sec) 

Yi  (sec-1) 

z1  (yd) 

1 

1.6468 

-0.01884 

0.0 

2 

1.6486 

0.418 

-95.897 

The  sound  speed  at  the  source  and  the  maximum  sound  speed  were,  respec- 
tively, 

c$  * 1.6471  kyd/sec 
^ v a 1.6486  kyd/sec 

maX 

resulting  in  a propagation  cone  of  2c  where 

c * Arccos  (c_/c  ) = 2.44°  . (4.5) 

5 maX 

Typically  we  ran  ACT.0R  using  e = 3.0°. 

After  making  several  calculations  with  ACTOR  using  the  two- 

layer  duct  we  noticed  a consistent  tendency  to  underestimate  the  Intensity 

at  the  50  foot  receiver.  A small  number  of  problems  were  subsequently 

(33 ) 

run  using  a three-layer  duct  suggested  by  Pedersen.'  ' This  was  the  form 
used  by  him  in  the  theoretical  calculations  reported  in  Reference  31. 

This  model  is  shown  in  Figure  4.2  with  the  relevant  parameters  being 
given  in  Table  4.3. 


Table  4.3.  Parameters  Characterizing  the  Three-Layer  Model. 


Layer 
Number,  i 

ci  (kyd/sec) 

Y.  (sec-1) 

zi  (ft) 

1 

1.64384 

0.062 

0 

2 

1.64322 

-0.015 

-30 

3 

1.64467 

0.390 

-321 

4.3  COMPARISON  OF  TRANSPORT  RESULTS  WITH  EXPERIMENT  AND  WITH 

THE  NORMAL  MOOE  RESULTS 

In  this  section  we  compare  the  transport  results  obtained  with 
the  ACTOR  code  with  the  experimental  results  and  with  those  obtained 
using  the  NOSC  normal  mode  code.  The  normal  mode  results  are  presented 
as  obtained  from  both  the  phased  addition  option  as  well  as  the  random 
phase  addition  option.  The  normal  mode  results,  as  well  as  the  experi- 
mental data,  were,  as  mentioned  in  the  previous  section,  provided  by 
M.  Pedersen  and  D.  White. 

The  transport  results  depend  upon  the  influence  of  two 
physical  mechanisms,  namely,  rough  surface  scattering  and  internal  wave 
scattering.  The  scattering  models  used  to  describe  these  mechanisms 
were  discussed  in  Section  3.  The  only  surface  scattering  model  used  was 
the  sawtooth  model  and  the  scattering  by  internal  waves  was  described 
by  the  Garrett-Munk  model. 

Before  presenting  any  of  the  results,  an  explanation  con- 
cerning the  captions  is  In  order.  The  slope  of  the  sawtooth  for  the 
surface  model  is  specified  by  the  symbol  S.  When  there  is  only 
specular  reflection  from  a flat  surface,  S * 0. 

The  internal  wave  scattering  is  characterized  by  one  of  the 
following  words  in  parentheses:  general,  linear,  or  quadratic.  The 
explanation  is  aided  by  considering  Figure  4.3.  Above  the  bottom  of 
the  surface  duct  three  options  are  shown  for  specifying  aJW  the 
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scattering  cross  section  due  to  internal  waves.  The  general  expression 
for  OjW  is  given  by  Eqs.  (3.91)  and  (3.92). 

The  first  calculation  run  with  internal  waves  used  the  general 
expression  for  The  results  are  shown  in  Figures  4.4  - 4.6  for 

v * 1030  Hz  at  the  receiver  depths  of  50,  232,  and  400  feet,  respectively. 
We  see  immediately  from  Figure  4.4  that  the  scattering  near  the  surface 
is  much  too  strong  since  the  transport  results  fall  progressi vely  farther 
below  the  experimental  data  as  range  increases.  This  is  not  surprising 
because  the  general  expression  does  not  apply  to  the  almost  isothermal 
conditions  which  exist  in  the  surface  duct.  In  fact,  we  expect  the 
scattering  cross  section  to  approach  zero  at  the  surface. 

We  picked  two  options  for  letting  aIW  go  to  zero  at  the 
surface.  They  are  shown  in  Figure  4.3  along  with  the  general  expression. 
The  linear  option  is  a simple  linear  interpolation  in  depth  between  zero 
at  the  surface  and  the  general  expression  value  at  the  bottom  of  the 
duct,  Zd>  The  quadratic  option  does  quadratic  interpolation  such  that 
OjW  = 0 at  the  surface,  aJW  * <jjW  (general  expression ) at  Z * Zd  and 
°IW  * ^ °IW  (^eral  expression)  at  Z * 3/4  Zd. 

The  case  where  S * 0 (specular  reflection  from  a flat  surface) 
and  with  the  linear  option  is  shown  in  Figures  4.7  - 4.9.  We  now  see 
that  the  transport  results  at  the  50  foot  receiver  are  in  much  better 
agreement  with  the  experimental  data.  Also,  the  transport  results  for 
the  232  foot  and  400  foot  receivers  are  quite  good.  At  short  range 
(<  10  kyd)  for  the  400  foot  cross  layer  case  there  is  some  discrepancy 
but,  with  the  diffraction  contribution,  the  calculated  and  experimental 
results  would  be  in  reasonably  good  agreement.  An  estimate  of  the 
diffraction  contribution  is  obtained  from  the  normal  mode  results.  At 
the  232  foot  receiver  the  transport  results  are  in  quite  good  agreement 
with  the  data.  To  be  sure,  the  transport  method  does  not  reproduce  all 
of  the  structure  observed  in  the  experimental  data.  This  is  primarily 
because  the  transport  equation  does  not  include  diffraction  effects  due 
to  the  deterministic  part  of  the  sound  speed  profile. 
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BILINERR  PROFILE,  S-O,  INTERNAL  WAVE  SCATTERING (GENERAL) 
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FIG.  4.6  , BILINEAR  PROFILE,  S-O,  INTERNAL  WAVE  SCATTERING (GENERAL) 


BILINEAR  PROFILE, S-O,  INTERNAL  NAVE  SCATTERING (LINEAR 
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FIG.  4.9  , BILINEAR  PROFILE, S-0,  INTERNAL  WAVE  SCATTERING (LINEAR! 


We  next  checked  to  see  how  the  transport  results  would  be 
modified  if  we  switched  to  the  quadratic  option,  still  using  S = 0.. 

The  results  are  shown  in  Figures  4.10  - 4.12.  The  transport  results 
for  the  50  foot  receiver  are  raised  considerably  and  now  agree  almost 
exactly  with  the  normal  mode  results.  The  transport  results  for  the 
232  foot  receiver  (Figure  4.11)  are  in  excellent  agreement  with  the 
experimental  data..  Unfortunately,  at  the  400  foot  depth  we  see  from 
Figure  4.12  that  not  enough  energy  is  being  scattered  out  of  the  duct 
to  this  below  duct  receiver.  This  is  rather  disappointing  in  view  of 
the  excellent  agreement  with  the  experimental  data  at  the  other  two 
depths. 

The  fact  that  the  structure  at  the  232  foot  level  was  predicted 
so  well  in  spite  of  the  neglect  of  diffraction  suggested  that  we  should 
increase  the  strength  of  the  internal  wave  scattering  by  a factor  of  two. 
In  the  duct  and  near  the  bottom  the  scattering  model  could  easily 
be  wrong  by  such  an  amount.  The  results,  however,  as  shown  in 
Figures  4.13  - 4.15  indicate  a relatively  weak  dependence  on  the  magnitude 
of  cr  when  the  quadratic  option  is  used.  We  note  a slight  deterioration 
of  the  transport  results  at  232  feet  and  practically  no  improvement  at 
the  400  foot  depth.  The  results  indicate  that  a substantial  change  in 
the  scattering  cross  section  would  have  to  be  made  within  the  framework 
of  the  quadratic  option  in  order  to  get  the  correct  amount  of  scattering 
below  the  duct.  We  would  then  expect  also  a further  deterioration  of 
the  232  foot  results. 

With  these  results  in  mind  we  next  turned  to  the  case  where 
only  surface  scattering  was  used.  Several  calculations  were  made  but 
only  the  results  of  the  best  one  will  be  presented.  These  are  contained 
in  Figures  4.16  - 4.18.  A similarity  is  observed  between  these  results 
and  those  shown  in  Figures  4.7  - 4.9  for  the  case  of  the  linear  option 
with  a flat  surface.  At  the  400  foot  level  (Figure  4.18)  we  see  what 
appear  to  be  shadow  zones  in  the  transport  results  at  about  7 and  14 
kyd.  This  prompted  us  to  combine  surface  scattering  (S  = 0.0025) 
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FIG.  4.11,  BILINEAR  PROFILE,  S-O,  INTERNAL  WAVE  SCATTER ING (QUADRAT IC1 
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FIG.  4.13,  BILINEAR  PROFILE, S-O, INTERNAL  WAVE  SCATTERING (2  X QUADRATIC) 
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FIG.  4.14,  BILINERR  PROFILE, S-0, INTERNAL  WAVE  SCATTERING (2  X QUADRATIC! 
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FIG.  4.15,  BILINEAR  PROFILE, S-0, INTERNAL  WAVE  SCATTERING (2  X QUADRAT I Cl 
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FIG.  4.16,  BILINEAR  PROFILE,  S-.0025,  NO  INTERNAL  WAVE  SCATTERING 
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first  with  the  linear  internal  wave  model  and  then  with  the  quadratic 
model.  These  results  are  given  in  Figures  4.19  - 4.24.  Observation  of 
Figure  4.21  indicates  that  the  linear  model  with  surface  scattering  Is  a 
little  too  strong.  The  quadratic  model  with  surface  scattering,  however, 
is  about  right  at  all  levels.  It  might  be  considered  somewhat  low  at  the 
50  foot  level,  but  this  level  Is  such  that  mode  interference  effects  could 
be  playing  a role.  Thus,  less  emphasis  is  placed  on  achieving  accurate 
agreement  with  the  experimental  data  here  than  at  the  two  lower  levels. 

The  results  for  the  S * 0.0025  quadratic  model  indicated 
that,  at  least  as  a sensitivity  check,  we  should  rerun  it  with  S * 

0.00125.  This  would  allow  us  to  examine  the  effect  of  a factor  of  two 
in  the  sawtooth  slope.  These  results  are  presented  in  Figures  4.25  - 4.27. 
We  immediately  see  that  the  energy  scattered  to  the  400  foot  level  Is  too 
small  by  a substantial  amount.  Hence,  we  conclude  that  the  below  duct 
ensonification  Is  much  more  sensitive  to  the  surface  scattering  parameter 
than  to  the  Internal  wave  parameter,  at  least  for  this  particular  problem. 

We  are  thus  left  with  some  uncertainty  in  how  to  ascribe  the 
below  duct  ensonification  to  the  various  scattering  mechanisms.  The 
situation  could  have  been  more  easily  clarified  if  accurate  surface 
condition  data  had  been  taken  and  if  better  range-depth  and  time- 
dependent  temperature  and  salinity  data  were  available.  It  is  clear, 
however,  that  either  a very  small  surface  swell  component  or  a reasonable 
internal  wave  model  In  the  duct  can  account  for  the  level  on  ensonification 
measured  at  the  below- layer  receiver. 

Before  going  on  to  consider  the  530  Hz  data,  we  want  to  look  at 
the  statistical  fluctuations  in  the  Monte  Carlo  data  due  to  using  a 
finite  number  of  rays  or  particles.  We  do  this  by  comparing  the  results 
in  Figures  4.28  - 4.30  with  those  in  Figures  4.22  - 4.24.  As  indicated 
in  the  captions,  the  calculation  used  to  generate  Figures  4.28  - 4.30 
used  7000  particles,  whereas  the  calculations  used  to  generate  all  the 
other  figures  used  either  3000  or  4000  particles.  We  see  that  the 
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FIG.  4.19,  BILINEAR  PROFILE,  S-.0025,  INTERNAL  WAVE  SCATTERING (LINEAR) 


FIG.  4.20,  BILINEAR  PROFILE,  S-.0025,  INTERNAL  WAVE  SCATTERING (LINEAR) 


FIG.  4.21,  BILINEAR  PROFILE,  S-.0025,  INTERNAL  WAVE  SCATTERING f LINEAR) 
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FIG.  4.22,  BILINEAR  PROFILE,  S-. 0025, INTERNAL  WAVE  SCATTERING (QUADRAT IC) 
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FIG.  4.30,  BILINEAR  PROFILE, S-. 0025,  1 NT.  WAVE  SCATTER  I NG I QUAO 1 7000  RAYS 


RflNeendLornRDS) 


smoothness  in  the  results  obtained  using  7000  particles  is  only  slightly 
better  than  when  3000  - 4000  particles  were  used.  After  the  variance 
calculation  technique  (discussed  in  Section  2.3.3)  is  incorporated  into 
ACTOR,  we  will  be  able  to  place  error  bars  on  the  curves. 

The  530  Hz  results,  both  experimental  and  calculated,  are  shown 
in  Figures  4.31  - 4.42.  At  530  Hz  scattering  by  internal  waves  is  much 
weaker  (by  a factor  of  four)  than  at  1030  Hz,  so  essentially  the  530  Hz 
data  will  only  be  used  as  further  confirmation  that  the  transport  results 
using  surface  scattering  (S  = .0025)  with  the  quadratic  internal  wave 
model  in  the  duct  are  consistent  with  the  experimental  data  and,  in  fact, 
provide  the  best  fit  at  530  Hz  as  well  as  at  1030  Hz. 

Starting  with  Figures  4.31  - 4.33  we  have  the  results  using 
only  rough  surface  scattering.  The  transport  results  are  in  reasonably 
good  agreement  with  the  data,  with  the  exception  of  the  short  range  400 
foot  case.  In  Figures  4.34  - 4.36  we  consider  a flat  surface  and  the 
linear  model  for  internal  wave  scattering  in  the  duct.  We  see  that 
without  some  rough  surface  scattering  the  intensity  is  overestimated 
at  the  50  foot  receiver.  The  232  foot*  results  are  not  bad,  and  the  short 
range  predictions  at  400  feet  are  better  than  with  rough  surface 
scattering  only.  This  suggests,  as  it  did  at  1030  Hz,  that  a combination 
of  rough  surface  scattering  and  internal  wave  scattering  be  used.  These 
results  are  displayed  in  Figures  4.37  - 4.42. 

The  normal  mode  calculations  do  not  include  either  rough 
surface  effects  or  internal  wave  effects.  Nevertheless,  we  see  from 
Figure  4.39  that  at  400  feet  the  normal  mode  results  are  in  good  agree- 
ment with  experiment,  indicating  that  at  530  Hz  diffraction  alone  can 
account  for  essentially  all  of  the  below-duct  ensonifi cation.  Thus,  as 
with  the  1030  Hz  case  we  again  select,  as  the  best  fit  to  the  experi- 
mental data,  the  results  obtained  using  a surface  roughness  parameter  of 
S 3 0.0025  and  the  quadratic  internal  wave  model  in  the  duct. 
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FIG.  4.33,  BILINEAR  PROFILE,  S-.0025,  NO  INTERNAL  HAVE  SCATTERING 
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FIG.  4.35,  BILINEAR  PROFILE,  S-0,  INTERNAL  WAVE  SCATTERING (LI NEAR) 
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FIG.  4.36,  BILINEAR  PROFILE,  S-O,  INTERNAL  WAVE  SCATTERING (LI NEAR) 


37,  BILINEAR  PROFILE,  S-.0025,  INTERNAL  WAVE  SCATTER ING I L I NEAR ) 
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FIG.  4.38,  BILINEAR  PROFILE,  S-.0025,  INTERNAL  WAVE  SCATTERING (LINEAR1 


40,  BILINEAR  PROFILE  , S«. 0025, INTERNAL  WAVE  SCATTER INGf QUADRATIC) 
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FIG.  4.41,  BILINEAR  PROFILE  , S-. 0025,  INTERNAL  WAVE  SCATTERING (QUADRATIC) 
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To  complete  the  transport  calculational  program,  we  calculated 
the  transport  solution  using  the  three-layer  model  discussed  in  Section 
4.2.  This  model  was  used  in  the  early  ray  calculations  described  in 
Reference  31.  We  decided  to  use  it  here  because  in  Figure  4.22  at  the 
50  foot  receiver  the  transport  results  appear  to  fall  below  the  experi- 
mental data.  As  mentioned  earlier  this  may  be  because  of  mode  interference 
effects  but  we  wanted  to  check  the  three-layer  model  for  model  sensitivity 
in  any  case.  As  can  be  seen  from  Figures  4.43  - 4.45  the  1030  Hz  transport 
results  at  the  50  foot  receiver  are  now  in  better  agreement  with  the 
experimental  data.  Also,  at  the  other  two  depths  agreement  with  experiment 
is  not  bad.  From  this  one  might  conclude  that  the  three-layer  model  is 
actually  a better  description  of  the  surface  duct. 

To  further  test  the  three-layer  model  we  also  ran  the  530  Hz 
case,  the  results  of  which  are  shown  in  Figures  4.46-4.48.  It  is 
immediately  clear  from  Figure  4.46  that  the  three-layer  model  predicts 
intensities  that  are  too  large  at  the  50  foot  receiver.  Because  the 
232  foot  results  are  somewhat  insensitive  to  whether  the  two-  or  three- 
layer  profile  is  used,  and  because  the  400  foot  results  are  probably 
dominated  by  diffraction,  we  can  only  say,  based  on  the  50  foot  data, 
that  the  two- layer  model  is  somewhat  better. 

We  had  to  state  that  the  400  foot  results  are  probably  domi- 
nated by  diffraction  because  we  actually  have  no  normal  mode  calcula- 
tions for  the  three-layer  model.  The  normal  mode  results  shown  in 
Figures  4.46  - 4.48  are  for  the  two-layer  model;  we  are  assuming  in 
making  the  statement  above  concerning  diffraction  that  it  will  be  rather 
insensitive  to  whether  the  two-  or  three-layer  model  is  used. 
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PIG.  4.43,  TRILINEAR  PROFILE, S-. 0025, INTERNAL  WAVE  SCATTER INGf QUADRATIC) 
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.45,  TRILINEAR  PROFILE,  S-.0025,  INTERNAL  HAVE  SCATTER I NG ( QUADRAT  I C 1 


.46,  TRILINE  FIR  PROFILE,  S-.  0025,  INTERNAL  WAVE  SCATTERING  I QUADRATIC) 
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FIG.  4.47,  TRILINEAR  PROFILE, S-. 0025, INTERNAL  WAVE  SCATTERING (QUADRATIC) 
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TRILINEAR  PROFILE, S-. 0025, INTERNAL  WAVE  SCATTERING (QUADRAT IC) 
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5.  CONCLUSIONS 


The  use  of  the  transport  equation  for  calculating  acoustic 
propagation  in  a medium  where  rough  surface  and  volume  scattering  effects 
are  Important  has  been  investigated.  The  results  of  the  study  Indicate 
that,  for  the  surface  duct  considered,  and  at  1030  Hz,  the  below-duct 
ensonifl cation  can  be  ascribed  to  two  effects.  The  first  of  these  Is 
scattering  by  a rough  surface.  No  detailed  surface  roughness  measure- 
ments were  made  during  the  acoustic  propagation  experiments  that  we-were^ — 
^comparing-eer, calculations*  with,  but  -eej* study  showed  that  even  a very 
small  swell  component  could  account  for  much  of  the  energy  scattered  below 
the  duct.  There  was  some  ambiguity, in^e/^ctly  how  to  describe  the  surface 
roughness.  -used,  a sawtooth  model  wWoJj  accounte^for  both  shadowing 
and  multiple  scatter.  Preliminary  comparisons  of  this  model  with  a more 
realistic  sinusoidal  model  suggest  that  little  difference  in  the  propa- 
gation results  would  have  been  observed  had  the  sinusoidal  model  actually 
been  used  in  the  calculations.^ 

The  Inclusion  of  volume  scattering  by  Internal  waves  was  also 
studied  using  the  Garrett-Munk  description.  It  was  observed  that  their 
prescription  for  the  scattering  cross  section  could  not  be  used  in  the 
surface  duct  Itself  because  It  led  to  a far  greater  scattering  of  energy 
out  of  the  duct  than  was  observed  experimentally.  The  more  physically 
reasonable  assumption  was  made  that  the  scattering  due  to  Internal  waves 
vanishes  at  the  surface.  Two  methods  of  approaching  zero  were  tried 
(linear  and  quadratic  Interpolation)  and  it  was  observed  that  the  propa- 
gation results  were  relatively  insensitive  to  exactly  how  the  cross 
section  approached  zero  at  the  surface  as  long  as  It  was  linear  or 
stronger. 

From  this  study  we  can  conclude  that  the  transport  equation  Is 
a powerful  tool  for  calculating  acoustic  propagation  In  a medium  where 
random  effects  are  Important.  It  is  clear,  however,  that  more  accurate 
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oceanographic  measurements  will  be  required  If  propagation  losses  are  to 
be  ascribed  to  a particular  scattering  mechanism.  A more  recent  experiment 
(SUDS  I)  has  been  completed  for  which  detailed  surface  roughness  and 
sound  speed  data  are  available.  The  application  of  the  transport  equation 
In  attempting  to  explain  this  data  would  be  a further  Interesting  test  of 
the  technique  and  also  of  value  to  the  underwater  sound  propagation 
community. 

Other  areas  where  the  transport  equation  could  fruitfully  be 
applied  are: 

1.  Anomalous  volume  attenuation  studies; 

2.  Convergence  zone  peaks/widths  considering  diffraction 
versus  scattering  in  spreading  convergence  zone 
boundaries; 

3.  Broadband  signal  propagation; 

4.  Frequency  spreading  in  a real  ocean  due  to  Doppler 
shifting;  and 

5.  Bottom  Interactions. 

In  short,  any  application  where  random  media  effects  are  Important  can  be 
efficiently  and  simply  addressed  using  the  transport  equation  technique. 
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APPENDIX  A 


PRELIMINARY  TEST  PROBLEM  FOR  THE  ACTOR  CODE 

Because  ACTOR  is  a totally  new  code  in  the  underwater  acoustic 
propagation  area,  we  felt  that  it  would  be  desirable  to  run  a test  problem 
without  surface  scattering  or  internal  wave  scattering  and  compare  the 
results  obtained  with  a standard  ray  theory  code  also  without  the  men- 
tioned scattering  mechanisms.  Under  these  conditions  the  two  codes  should 
be  in  almost  exact  agreement.  The  problem  selected  for  this  purpose  was 
the  two-layer  duct  whose  properties  were  given  in  Table  4.2. 

As  was  the  case  during  the  actual  experiment  discussed  in 
Section  4.1  the  source  was  at  a depth  of  55  feet  and  the  propagation  loss 
versus  range  curves  were  calculated  at  depths  of  50,  232  and  400  feet. 

A comparison  of  propagation  loss  as  calculated  by  the  ray  theory  code  and 
the  transport  code  for  these  depths  is  given  in  Figures  A.l  - A. 3.  The 
corresponding  ray  diagram  is  provided  in  Figure  A. 4.  Both  the  ray  diagram 
and  the  ray  code  propagation  loss  data  were  provided  by  D.  White  andf 
M.  Pedersen*  of  NOSC. 

These  figures  clearly  indicate  that  the  ray  and  transport  codes 
are  in  very  good  agreement.  Also,  from  the  ray  diagram,  we  see  that  at 
the  232  foot  level  the  transport  code  is  predicting  shadow  zones  at  the 
right  places  and  of  the  right  width  in  range. 

Also,  at  the  400  depth,  both  codes  predict  a total  shadow  zone 
beyond  5 kyd.  However,  as  mentioned  in  Section  4,  this  region  can  be 
ensonlfled  due  to  scattering  or  diffraction. 

The  calculations  discussed  in  this  Appendix  were  actually  run 
before  any  attempt  was  made  to  compare  calculated  results  with  experi- 
mental data.  The  results  of  this  test  problem  indicated  that  the  trans- 
port code  was  doing  the  particle  tracking  correctly  and  that  the  absorption 
and  editing  algorithms  were  properly  formulated. 


* White,  0.  and  Pedersen,  M.,  private  communication. 


Figure  A. 2.  Propagation  Loss  Versus  Range  for  the  232  ft  Receiver 
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Figure  A, 4.  Ray  Diagram  for  the  Two-Layer  Duct  Showing  the  Shadow 
Zones  for  the  232  Foot  Receiver.  The  Limiting  Rays 
for  Propagation  are  at  ±2.41°. 
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APPENDIX  B 

RANDOM  NUMBER  GENERATION  TEST 

The  random  number  generation  scheme  used  in  ACTOR  makes  use  of 
the  multiplicative  congruential  method.*  We  ran  singlet,  doublet,  and 
triplet  tests,  each  type  consisting  of  three  sets. 

The  random  number  generator  is  given  by 

48 

x = (mod  m) , m = 2 

The  set  2 x is  the  last  random  number  from  set  1 and  the  set  3 xrt  is 
O u 

the  last  random  number  from  set  2.  The  parameters  for  the  various  sets 
were  as  given  in  Table  B.l. 

Table  B.l.  Parameters  for  Uniformity  Tests. 

SET  1 SET  2 SET  3 

Xo=20001057152522201031  xo=20007050346512654371  xo=20006443541674732631 

X-20000000000 100000053  X=same  as  set  1 X=same  as  set  1 


; 


t 


The  chi-square  test  is  expressed  by  the  following  equation* 


A ivs  - "Ps)2 
■ t*l  np 


(B.l) 


+ Spanier,  J.  and  Gelsord,  E.  M. , Monte  Carlo  Principles  and  Neutron 
Transport  Problems,  Addison  Wesley  Publishers,  (1969). 

* Knuth,  D.  E.,  The  Art  of  Computer  Programming,  Volume  2,  Addison-Wesley , 
Reading  Massachusetts  (1969). 
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where 

n = number  of  independent  trials 

k = number  of  possible  outcomes  of  a single  trial 
Y$  = actual  number  of  trials  which  resulted  in  outcome  s 

Ps  = probability  that  each  trial  results  in  outcome  s. 

We  ran  the  random  number  test  so  that  there  were  729  possible  outcomes 
for  each  set  and  n was  set  equal  to  729,000.  The  singlet  test  places 
each  random  number  in  the  appropriate  equi-sized  slot  in  a one-dimensional 
array  of  size  729.  The  doublet  test  places  successive  pairs  of  random 
numbers  in  the  appropriate  equi-sized  slots  in  a 27  x 27  two-dimensional 
array  and  729,000  pairs  were  generated.  The  triplet  test  proceeds  simi- 
larly for  a three-dimensional  array  of  size  9x9x9.  To  fully  test  a 
random  number  generator  all  p-tuples  should  be  tested  where  p -*■  ». 

Clearly  this  is  a practical  impossibility  and  we  have  settled  on  p 5 3. 

The  results  of  all  these  tests  are  shown  in  Figure  B.l.  Values 
of  V from  Eq.  (B.l)  which  lie  above  820  or  below  640  would  be  suspect. 

This  would  especially  be  so  if  it  happened,  say,  two  out  of  three  times. 

As  can  be  seen  from  the  figure,  however,  for  all  of  our  tests  the  value 
of  V lies  in  the  acceptable  range.  In  this  figure  0^  refers  to  the 
first  singlet  set,  O2  the  second,  etc. 

We  have  also  included  the  results  obtained  using  the  Control 
Data  supplied  random  number  generator.  Only  one  set  of  each  type  was 
checked  for  the  CDC  case. 

As  can  be  seen  from  the  figure,  both  random  number  generators 
appear  to  be  satisfactory.  From  these  results  one  might  ask  why  we 
bothered  with  our  own  random  number  generator  since  the  one  provided  by 
CDC  appears  to  be  adequate.  The  reason  is  that  in  going  from  one  instal- 
lation to  another  it  Is  never  certain  that  the  starting  seeds  are  the 
same.  In  order  that  certain  test  problems  be  reproducible,  we  have  found 
it  much  easier  to  carry  our  own  random  number  generator. 
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